Package osgeo :: Package utils :: Module gdalmove
[hide private]
[frames] | no frames]

Source Code for Module osgeo.utils.gdalmove

  1  #!/usr/bin/env python3 
  2  # ****************************************************************************** 
  3  # 
  4  #  Project:  GDAL Python Interface 
  5  #  Purpose:  Application for "warping" an image by just updating it's SRS 
  6  #            and geotransform. 
  7  #  Author:   Frank Warmerdam, warmerdam@pobox.com 
  8  # 
  9  # ****************************************************************************** 
 10  #  Copyright (c) 2012, Frank Warmerdam 
 11  # 
 12  #  Permission is hereby granted, free of charge, to any person obtaining a 
 13  #  copy of this software and associated documentation files (the "Software"), 
 14  #  to deal in the Software without restriction, including without limitation 
 15  #  the rights to use, copy, modify, merge, publish, distribute, sublicense, 
 16  #  and/or sell copies of the Software, and to permit persons to whom the 
 17  #  Software is furnished to do so, subject to the following conditions: 
 18  # 
 19  #  The above copyright notice and this permission notice shall be included 
 20  #  in all copies or substantial portions of the Software. 
 21  # 
 22  #  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 
 23  #  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
 24  #  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
 25  #  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
 26  #  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
 27  #  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 
 28  #  DEALINGS IN THE SOFTWARE. 
 29  # ****************************************************************************** 
 30   
 31  import math 
 32  import sys 
 33   
 34  from osgeo import gdal, osr 
 35   
 36  ############################################################################### 
 37   
 38   
39 -def fmt_loc(srs_obj, loc):
40 if srs_obj.IsProjected(): 41 return '%12.3f %12.3f' % (loc[0], loc[1]) 42 return '%12.8f %12.8f' % (loc[0], loc[1])
43 44 ############################################################################### 45 46
47 -def move(filename, t_srs, s_srs=None, pixel_threshold=None):
48 49 # ------------------------------------------------------------------------- 50 # Open the file. 51 # ------------------------------------------------------------------------- 52 ds = gdal.Open(filename) 53 54 # ------------------------------------------------------------------------- 55 # Compute the current (s_srs) locations of the four corners and center 56 # of the image. 57 # ------------------------------------------------------------------------- 58 corners_names = [ 59 'Upper Left', 60 'Lower Left', 61 'Upper Right', 62 'Lower Right', 63 'Center'] 64 65 corners_pixel_line = [ 66 (0, 0, 0), 67 (0, ds.RasterYSize, 0), 68 (ds.RasterXSize, 0, 0), 69 (ds.RasterXSize, ds.RasterYSize, 0), 70 (ds.RasterXSize / 2.0, ds.RasterYSize / 2.0, 0.0)] 71 72 orig_gt = ds.GetGeoTransform() 73 74 corners_s_geo = [] 75 for item in corners_pixel_line: 76 corners_s_geo.append( 77 (orig_gt[0] + item[0] * orig_gt[1] + item[1] * orig_gt[2], 78 orig_gt[3] + item[0] * orig_gt[4] + item[1] * orig_gt[5], 79 item[2])) 80 81 # ------------------------------------------------------------------------- 82 # Prepare a transformation from source to destination srs. 83 # ------------------------------------------------------------------------- 84 if s_srs is None: 85 s_srs = ds.GetProjectionRef() 86 87 s_srs_obj = osr.SpatialReference() 88 s_srs_obj.SetFromUserInput(s_srs) 89 90 t_srs_obj = osr.SpatialReference() 91 t_srs_obj.SetFromUserInput(t_srs) 92 93 tr = osr.CoordinateTransformation(s_srs_obj, t_srs_obj) 94 95 # ------------------------------------------------------------------------- 96 # Transform the corners 97 # ------------------------------------------------------------------------- 98 99 corners_t_geo = tr.TransformPoints(corners_s_geo) 100 101 # ------------------------------------------------------------------------- 102 # Compute a new geotransform for the image in the target SRS. For now 103 # we just use the top left, top right, and bottom left to produce the 104 # geotransform. The result will be exact at these three points by 105 # definition, but if the underlying transformation is not affine it will 106 # be wrong at the center and bottom right. It would be better if we 107 # used all five points for a least squares fit but that is a bit beyond 108 # me for now. 109 # ------------------------------------------------------------------------- 110 ul = corners_t_geo[0] 111 ur = corners_t_geo[2] 112 ll = corners_t_geo[1] 113 114 new_gt = (ul[0], 115 (ur[0] - ul[0]) / ds.RasterXSize, 116 (ll[0] - ul[0]) / ds.RasterYSize, 117 ul[1], 118 (ur[1] - ul[1]) / ds.RasterXSize, 119 (ll[1] - ul[1]) / ds.RasterYSize) 120 121 inv_new_gt = gdal.InvGeoTransform(new_gt) 122 123 # ------------------------------------------------------------------------- 124 # Report results for the five locations. 125 # ------------------------------------------------------------------------- 126 127 corners_t_new_geo = [] 128 error_geo = [] 129 error_pixel_line = [] 130 corners_pixel_line_new = [] 131 132 print('___Corner___ ________Original________ _______Adjusted_________ ______ Err (geo) ______ _Err (pix)_') 133 134 for i in range(len(corners_s_geo)): # pylint: disable=consider-using-enumerate 135 136 item = corners_pixel_line[i] 137 corners_t_new_geo.append( 138 (new_gt[0] + item[0] * new_gt[1] + item[1] * new_gt[2], 139 new_gt[3] + item[0] * new_gt[4] + item[1] * new_gt[5], 140 item[2])) 141 142 error_geo.append((corners_t_new_geo[i][0] - corners_t_geo[i][0], 143 corners_t_new_geo[i][1] - corners_t_geo[i][1], 144 0.0)) 145 146 item = corners_t_geo[i] 147 corners_pixel_line_new.append( 148 (inv_new_gt[0] + item[0] * inv_new_gt[1] + item[1] * inv_new_gt[2], 149 inv_new_gt[3] + item[0] * inv_new_gt[4] + item[1] * inv_new_gt[5], 150 item[2])) 151 152 error_pixel_line.append( 153 (corners_pixel_line_new[i][0] - corners_pixel_line[i][0], 154 corners_pixel_line_new[i][1] - corners_pixel_line[i][1], 155 corners_pixel_line_new[i][2] - corners_pixel_line[i][2])) 156 157 print('%-11s %s %s %s %5.2f %5.2f' % 158 (corners_names[i], 159 fmt_loc(s_srs_obj, corners_s_geo[i]), 160 fmt_loc(t_srs_obj, corners_t_geo[i]), 161 fmt_loc(t_srs_obj, error_geo[i]), 162 error_pixel_line[i][0], 163 error_pixel_line[i][1])) 164 165 print('') 166 167 # ------------------------------------------------------------------------- 168 # Do we want to update the file? 169 # ------------------------------------------------------------------------- 170 max_error = 0 171 for err_item in error_pixel_line: 172 this_error = math.sqrt(err_item[0] * err_item[0] + err_item[1] * err_item[1]) 173 if this_error > max_error: 174 max_error = this_error 175 176 update = False 177 if pixel_threshold is not None: 178 if pixel_threshold > max_error: 179 update = True 180 181 # ------------------------------------------------------------------------- 182 # Apply the change coordinate system and geotransform. 183 # ------------------------------------------------------------------------- 184 if update: 185 ds = None 186 ds = gdal.Open(filename, gdal.GA_Update) 187 188 print('Updating file...') 189 ds.SetGeoTransform(new_gt) 190 ds.SetProjection(t_srs_obj.ExportToWkt()) 191 print('Done.') 192 193 elif pixel_threshold is None: 194 print('No error threshold in pixels selected with -et, file not updated.') 195 196 else: 197 print("""Maximum check point error is %.5f pixels which exceeds the 198 error threshold so the file has not been updated.""" % max_error) 199 200 ds = None
201 202 ############################################################################### 203 204
205 -def Usage():
206 print(""" 207 gdalmove.py [-s_srs <srs_defn>] -t_srs <srs_defn> 208 [-et <max_pixel_err>] target_file 209 """) 210 sys.exit(1)
211 212
213 -def main(argv):
214 # Default GDAL argument parsing. 215 216 argv = gdal.GeneralCmdLineProcessor(argv) 217 if argv is None: 218 sys.exit(0) 219 220 if len(argv) == 1: 221 Usage() 222 223 # Script argument defaults 224 s_srs = None 225 t_srs = None 226 filename = None 227 pixel_threshold = None 228 229 # Script argument parsing. 230 231 i = 1 232 while i < len(argv): 233 234 if argv[i] == '-s_srs' and i < len(argv) - 1: 235 s_srs = argv[i + 1] 236 i += 1 237 238 elif argv[i] == '-t_srs' and i < len(argv) - 1: 239 t_srs = argv[i + 1] 240 i += 1 241 242 elif argv[i] == '-et' and i < len(argv) - 1: 243 pixel_threshold = float(argv[i + 1]) 244 i += 1 245 246 elif filename is None: 247 filename = argv[i] 248 249 else: 250 print('Urecognised argument: ' + argv[i]) 251 Usage() 252 253 i = i + 1 254 # next argument 255 256 if filename is None: 257 print('Missing name of file to operate on, but required.') 258 Usage() 259 260 if t_srs is None: 261 print('Target SRS (-t_srs) missing, but required.') 262 Usage() 263 264 move(filename, t_srs, s_srs, pixel_threshold)
265 266 267 if __name__ == '__main__': 268 sys.exit(main(sys.argv)) 269