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())--若输出为空字符串则无投影。