python投影Himawari-8数据

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的影像如下: