|||
题记: 写该博客是因为自己经常遇到这个问题,而我发现网络上关于这方面浏览量高的一些代码竟然都有误,每次照搬都被虐得很惨。有一些同志在某些博客下方留言说代码有问题,而博主没有回应,也没有更改错误。为了自己及他人以后的方便,故于此附上我反复验证后的代码,并分享一下我解决这个问题时的一些经验。
摘要:(1)图像坐标——投影坐标——经纬度坐标 的映射代码
(2)经纬度坐标——投影坐标——图像坐标 的映射代码
(3)运行出现 NotImplementedError: Wrong number or type of arguments for overloaded
function 'CoordinateTransformation_TransformPoint'. 错误的解决办法。
正文:(1)图像坐标——投影坐标——经纬度坐标
Tiles_file='S2A_MSIL1C_20170309T031621_N0204_R075_T50SLG_20170309T031619.tif'
dataset = gdal.Open(Tiles_file) # 读入影像文件
nXSize,nYSize=dataset.RasterXSize,dataset.RasterYSize # 影像列数,行数
trans = dataset.GetGeoTransform()
prosrs = osr.SpatialReference()
prosrs.ImportFromWkt(dataset.GetProjection())
geosrs = prosrs.CloneGeogCS()
col=0 # 【图像坐标】像素所在列
row=0 # 【图像坐标】像素所在行
px = trans[0] + col * trans[1] + row * trans[2] # 299980.0 该点在经度方向的【投影坐标】
py = trans[3] + col * trans[4] + row * trans[5] # 4090220.0 该点在纬度方向的【投影坐标】
ct = osr.CoordinateTransformation(prosrs, geosrs) # 表示从投影坐标系映射到地理坐标系
coords = ct.TransformPoint(px, py)
lat=coords [0] # 【经纬度坐标】纬度
lon=coords [1] # 【经纬度坐标】经度
补充说明:
1. 如果你想要得到左上角像素的中心坐标,你需要设置 col,row=0.5,0.5
下面附上该影像的四至点坐标:
col,row=0,0 #299980.0 4200010.0 左上 像素左上角点 # (37.92566740942532, 114.7242415990631)
col,row=nXSize-1,nYSize-1 #409790.0 4090220.0 右下像素左上角点 # (36.953738446008366, 115.98676016473146)
col,row=nXSize-1,0 #409790.0 4200010.0 右上像素左上角点 # (37.94320029660642, 115.97332489462718)
col,row=0,nYSize-1 #299980.0 4090220.0 左下像素左上角点 # (36.936817812938365, 114.75399815989871)
2. 关于trans = dataset.GetGeoTransform() 的介绍
trans 是一个Tuple类型的变量包含6个元素,[trans[0],trans[3]] 是影像左上角像素的左上角点坐标,前者为经度方向的投影坐标,后者为纬度方向的投影坐标 trans[1] 表示像素宽度 ,trans[5] 表示像素高度
3. 关于coords = ct.TransformPoint(px, py) 的介绍
coords 是一个Tuple类型的变量包含3个元素,coords [0]为纬度,coords [1]为经度,coords [2]为高度
(2)经纬度坐标——投影坐标——图像坐标
ct2 = osr.CoordinateTransformation(geosrs, prosrs) # 表示从地理坐标系投映射到影坐标系
coords = ct2.TransformPoint(36.936817812938365, 114.75399815989871) # 【经纬度坐标】
Pro_lon=coords [0] # 【投影坐标】经度方向
Pro_lat=coords [1] # 【投影坐标】纬度方向
trans = dataset.GetGeoTransform()
col_id = int((Pro_lon - trans[0]) / trans[1]+0.000000001) # 【图像坐标】列方向 +0.000000001为了填补浮点数计算误差
row_id = int((trans[3] - Pro_lat ) / -trans[5]+0.000000001) #【图像坐标】行方向 -trans[5] 为正数
(3) 出现如下错误的解决方法
NotImplementedError: Wrong number or type of arguments for overloaded function'CoordinateTransformation_TransformPoint'. Possible C/C++ prototypes are: OSRCoordinateTransformationShadow:TransformPoint(double [3]) OSRCoordinateTransformationShadow::TransformPoint(double [4]) OSRCoordinateTransformationShadow::TransformPoint(double [3],double,double,double) OSRCoordinateTransformationShadow::TransformPoint(double [4],double,double,double,double)
https://gis.stackexchange.com/questions/331908/notimplementederror-wrong-number-or-type-of-arguments-for-overloaded-function建议先检查projection 是否为空,我print(dataset.GetProjection())
输出的果然是空字符串' '。
我最终的解决办法是:首先卸载gdal ,pip uninstall gdal(site-pakage 下若有同名残留文件,且删除)
然后,去下载gdal whl 文件 https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
GDAL‑3.0.4‑cp36‑cp36m‑win_amd64.whl 表示指出python3.6, window 64位系统
该文件下载很慢,如有需要发QQ邮件给我552569054@qq.com 最后,pip install GDAL‑3.0.4‑cp36‑cp36m‑win_amd64.whl
结语:如有错误,欢迎大家在下方留言并发邮件提醒我,谢谢。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-5-19 23:57
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社