popopochild的个人博客分享 http://blog.sciencenet.cn/u/popopochild

博文

Python GDAL 图像坐标,投影坐标,经纬度坐标 三者映射及运行错误解决

已有 9352 次阅读 2020-3-14 18:04 |个人分类:编程|系统分类:科研笔记| Python, GDAL, 图像坐标, 投影坐标, 经纬度坐标

题记: 写该博客是因为自己经常遇到这个问题,而我发现网络上关于这方面浏览量高的一些代码竟然都有误,每次照搬都被虐得很惨。有一些同志在某些博客下方留言说代码有问题,而博主没有回应,也没有更改错误。为了自己及他人以后的方便,故于此附上我反复验证后的代码,并分享一下我解决这个问题时的一些经验。


摘要:(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 


结语:如有错误,欢迎大家在下方留言并发邮件提醒我,谢谢。

                                        




https://wap.sciencenet.cn/blog-2991632-1223505.html

上一篇:MODIS LAI 产品关于时间的命名规则及合成法则
收藏 IP: 223.68.48.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (4 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-5-19 23:57

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部