1 import rasterio
2 from pyproj import Proj, transform, CRS
3
4 ### geotiff_file: the original remote sensing file contain latitude and longitude information
5 def coord_init(geotiff_file):
6 with rasterio.open(geotiff_file) as r:
7 # T0 = r.transform # upper-left pixel corner affine transform
8 up_left_east, up_left_north = r.transform * (0, 0)
9
10 bottom_right_east, bottom_right_north = r.transform * (r.width, r.height)
11 print('(up_left_east, up_left_north)', (up_left_east, up_left_north))
12 print('(bottom_right_east, bottom_right_north)', (bottom_right_east, bottom_right_north))
13
14 res = r.res# 分辨率
15
16 p0 = Proj(r.crs)
17 p1 = Proj(CRS.from_epsg(4326))
18 p2 = Proj(CRS.from_epsg(3857))
19
20 up_left_lat, up_left_lon = transform(p0, p1, up_left_east, up_left_north)# p0 --> p1
21 # x_east, y_north = transform(p1, p0, up_left_lat, up_left_lon)# p1 --> p0
22 # bottom_right_lat, bottom_right_lon = transform(p0, p1, bottom_right_east, bottom_right_north) # p0 --> p1
23 # x_east, y_north = transform(p1, p2, 39.258, 120.445)# p1 --> p2
24 # y_lat, x_lon = transform(p2, p1, x_east, y_north) # p2 --> p1
25 # x_botm_east, y_botm_north = transform(p1, p0, bottom_right_lat, bottom_right_lon) # p1 --> p0
26
27 return p0, p1, p2, (up_left_lon, up_left_lat), res
28
29 # input: points: [0]-> latitude/y, [1]-> longitude/x
30 # output: latitudes, longitudes
31 def coord_trans_array(p0, p1, p2, points, up_left_lonlat, res):
32 up_left_lon, up_left_lat = up_left_lonlat[0], up_left_lonlat[1]
33 X, Y = transform(p1, p2, up_left_lat, up_left_lon)# p1 --> p2
34
35 pix_x, pix_y = res[0], res[1]
36 eastings = []# x, longitude
37 northings = []# y, latitude
38 for p in points:
39 p_x = X + p[1] * abs(pix_x)
40 p_y = Y - p[0] * abs(pix_y)
41 eastings.append(p_x)
42 northings.append(p_y)
43 lats, longs = transform(p2, p1, eastings, northings)# p2 --> p1
44 x_east, y_north = transform(p1, p0, lats, longs) # p1 --> p0
45 points_new = list(zip(lats, longs))
46 return points_new
47
48 # input: point: [0]-> latitude/y, [1]-> longitude/x
49 # output: latitudes, longitudes
50 def coord_trans_point(p0, p1, p2, point, up_left_lonlat, res):
51 # point[0] -->north, point[1] -->east
52 # p0: input primitive EPSG
53 # p1: EPSG:4326
54 # p2: EPSG:3857
55 up_left_lon, up_left_lat = up_left_lonlat[0], up_left_lonlat[1]
56 X, Y = transform(p1, p2, up_left_lat, up_left_lon)# p1 --> p2
57 pix_x, pix_y = res[0], res[1]
58 X += point[1]*abs(pix_x)
59 Y -= point[0]*abs(pix_y)
60
61 lats, longs = transform(p2, p1, X, Y)# p2 --> p1
62 x_east, y_north = transform(p1, p0, lats, longs) # p1 --> p0
63 return [lats, longs]
64
65 if __name__ == '__main__':
66 geotiff_file = '/home/jiangshan/Projects/data/tif_input/10378780_15.tiff'
67 p0, p1, p2, up_left_lonlat, res = coord_init(geotiff_file)
68 points = [(0, 0), (0, 750), (0, 1500)]# y, x
69 pts = coord_trans_array(p0, p1, p2, points, up_left_lonlat, res)
70 print(pts)
71 pts = coord_trans_point(p0, p1, p2, (0, 750), up_left_lonlat, res)
72 print(pts) 个人学习记录