IDL是可以构建GLT对影像进行校正的,在python里面应该是用gdal.warp()实现的,据传可以编写VRT文件写入经纬度文件的路径就可以实现和GLT校正同样的功能,可惜我太菜了不会搞,就弄个笨办法来用了。
输入文件为HSD数据提取的影像或者NC文件中得到的影像数据,具体长这样:
import gdal
inputfile='*\\HS_H08_20190916_0300_FLDK_B3_R20.dat'
outtif='*********.tif'
memDs = gdal.Open(inputfile)
cols = memDs.RasterXSize
rows = memDs.RasterYSize
if rows == 11000 and cols == 11000:
res = 0.01
else:
res = 0.02
# 几何校正
# 定义空间参考
srs = osr.SpatialReference()
# 定义地球长半轴a=6378137.0m,地球短半轴b=6356752.3m,卫星星下点所在经度140.7,目标空间参考
srs.ImportFromProj4('+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=140.7 +no_defs')
memDs.SetProjection(srs.ExportToWkt())
memDs.SetGeoTransform([-cols*res*50000, int(res*100000), 0, cols*res*50000, 0, int(-res*100000)])
dstFilePath = os.path.join(outfolder,'H8_'+str(cols)+'.tif')
if os.path.exists(outtif):
os.remove(outtif)
warpDs = gdal.Warp(outtif, memDs, dstSRS='EPSG:4326', outputBounds=(60.0, -90.0, 222.0, 90.0), xRes=res, yRes=res)
del warpDs
(60.0, -90.0, 222.0, 90.0)的范围是利用IDL计算影像像元坐标时算出来的经纬度范围估算的外边界范围,用起来误差还是不太大的。校正完影像长这样(跟矢量很贴合了):
放大看这里:
tips:
srs.ImportFromProj4(’+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=140.7 +no_defs’)
代码中的数据换成FY4A的就可以对FY4A影像进行投影的,亲测可用。
outputBounds=(60.0, -90.0, 222.0, 90.0)替换为outputBounds=(80.0, -60.0, 200.0, 60.0),生成的影像即为NC文件对应的经纬度范围内的影像(省去裁剪了),不过NC影像大小为6001X6001,使用gdal.Warp生成的为6000X6000,问题不大,6000X6000的影像如下: