代码收藏家技术教程 2024-06-08

Python GDAL：TIF图像坐标、投影坐标和经纬度转换（解析i和j的顺序）

【太长不看版】

tif是一张具有地理编码信息的图像。

gdal能读出里面的基本地理编码信息，根据这6个基本的信息，可以推导出其他需要的地理信息。

trans(0) 左上像素左上角的x坐标。
trans(1) w-e像素分辨率/像素宽度。
trans(2) 行旋转（通常为零）。
trans(3) 左上像素左上角的y坐标。
trans(4) 列旋转（通常为零）。
trans(5) n-s像素分辨率/像素高度（北上图像为负值）

【详细举例版】

nXSize = dataset.RasterXSize #列数
nYSize = dataset.RasterYSize #行数

img_lons = np.zeros((nYSize, nXSize))
img_lats = np.zeros((nYSize, nXSize))

trans = dataset.GetGeoTransform()
prosrs = osr.SpatialReference()
prosrs.ImportFromWkt(dataset.GetProjection())
geosrs = prosrs.CloneGeogCS()。
ct = osr.CoordinateTransformation(prosrs, geosrs)


for col in range(nXSize):   
    for row in range(nYSize):
        px = trans[0] + col * trans[1] + row * trans[2]    # lons
        py = trans[3] + col * trans[4] + row * trans[5]    # lats
        
        
        coords = ct.TransformPoint(px, py)
        
        lat=coords[0]        
        lon=coords[1]
        
        img_lons[row,col] =  lon           
        img_lats[row,col] =  lat

网上很多代码都是错的，将col和row的位置弄反了。

实际上并不复杂。

以这个公式为例“px = trans[0] + col * trans[1] + row * trans[2]”，

trans[2]是旋转角度，通常为0，因此公式变为“px = trans[0] + col * trans[1]”，

px是我们要计算的x轴方向上的坐标，而col代表列数（XSize），trans[1]是x方向的分辨率。

换句话说，col是x方向上的间隔个数，而trans[1]是x方向上每个间隔的宽度。

因此 “px = trans[0] + col * trans[1]” → 每个px=原点+间隔数*间隔宽度。

（如果我讲清楚了，可以给我一个小赞赞吗，谢谢哈哈哈(∩_∩)）

【注意tif本身是否带有投影信息】

上一段代码先将图像坐标（即tran的6个参数）→转化为了投影坐标（即px、py）再→转为经纬度（lon、lat）

这是对于有投影的tif图像的处理；若tif本身没有投影信息，则以上代码算出的px、py即为经纬度。

检查投影信息方法：print（dataset.GetProjection()）–若输出为空字符串则无投影。