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