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)