1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31 import math
32 import sys
33
34 from osgeo import gdal, osr
35
36
37
38
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
51
52 ds = gdal.Open(filename)
53
54
55
56
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
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
97
98
99 corners_t_geo = tr.TransformPoints(corners_s_geo)
100
101
102
103
104
105
106
107
108
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
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)):
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
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
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
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
214
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
224 s_srs = None
225 t_srs = None
226 filename = None
227 pixel_threshold = None
228
229
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
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