tiff遥感图像空间坐标转换（工作太忙，仅仅作为记录） 

 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)

 

posted @ 2021-09-20 20:00  土博姜山山  阅读(610)  评论(0)  编辑  收藏  举报