风云三号极轨卫星的几何校正
风云三号极轨卫星的几何校正
ytkz遇到问题记录：
1.风云极轨卫星的几何校正，有单独的经度和纬度以及DN数据，试过GLT校正但是不行
解决方法：
构造vrt文件，再进行几何校正。
vrt文件
官方vrt资料
VRT驱动程序是GDAL的一个格式驱动程序，它允许虚拟GDAL数据集由其他GDAL数据集组成，并具有重新定位、可能应用的算法以及更改或添加的各种元数据。数据集的VRT描述可以保存为XML格式，通常给出扩展名.VRT。
针对风云三号提取EV_Emissive、EV_RefSB数据示例
```python
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2021/12/24 10:21
# @Author :
# @File : fy3_corretion.py
from osgeo import gdal, osr
import os
def geoMERSI2(file, geoFile, afterGeoPath):
'''
Parameters
----------
file : 文件绝对路径
需要校正的文件.
geoFile : 文件绝对路径
地理坐标存在的文件.
afterGeoPath : tif文件夹绝对路径
经过校正后的tif格式 Albers投影 250m空间分辨率.
Returns
-------
None.
'''
# os.chdir(r'F:\Py_project\gdalDealMERSI2')
dataset = gdal.Open(file)
# geoDataset = gdal.Open(geoFile)
# print(gdal.Info(dataset)) # 查看元数据
# print(gdal.Info(geoDataset))
# subdataset = dataset.GetSubDatasets() # 查看子数据集
# print(subdataset)
'''
查看子数据集位置：
SUBDATASET_6_NAME=HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_0250M_MS.hdf"://Data/EV_250_Emissive_b24
SUBDATASET_6_DESC=[8000x8192] //Data/EV_250_Emissive_b24 (16-bit unsigned integer)
'''
'''
地理校正数据的位置：
SUBDATASET_1_NAME=HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Latitude
SUBDATASET_1_DESC=[8000x8192] //Latitude (32-bit floating-point)
SUBDATASET_2_NAME=HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Longitude
SUBDATASET_2_DESC=[8000x8192] //Longitude (32-bit floating-point)
'''
vrtEV_EmissiveDir = os.path.splitext(file)[0]+'EV_Emissive' + '.vrt'
vrtEV_RefSBDir = os.path.splitext(file)[0] + 'EV_RefSB' + '.vrt'
subDatasetEV_Emissive = dataset.GetSubDatasets()[0][0] #
subDatasetEV_RefSB= dataset.GetSubDatasets()[1][0]
vrtEV_EmissiveFile = gdal.Translate(vrtEV_EmissiveDir,
subDatasetEV_Emissive,
format='vrt')
vrtEV_RefSBDirFile = gdal.Translate(vrtEV_RefSBDir,
subDatasetEV_RefSB,
format='vrt')
'''
需要写入描述GEOLOCATION元数据域的信息：
<Metadata domain="GEOLOCATION">
<MDI key="LINE_OFFSET">1</MDI>
<MDI key="LINE_STEP">1</MDI>
<MDI key="PIXEL_OFFSET">1</MDI>
<MDI key="PIXEL_STEP">1</MDI>
<MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</MDI>
<MDI key="X_BAND">1</MDI>
<MDI key="X_DATASET">HDF5:"D:\dengkaiyuan\data\FY3C_VIRRX_GBAL_L1_20211110_0825_GEOXX_MS.HDF"://Geolocation/Longitude</MDI>
<MDI key="Y_BAND">1</MDI>
<MDI key="Y_DATASET">HDF5:"D:\dengkaiyuan\data\FY3C_VIRRX_GBAL_L1_20211110_0825_GEOXX_MS.HDF"://Geolocation/Latitude</MDI>
</Metadata>
'''
srs = osr.SpatialReference()
srs.ImportFromProj4(
'+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
lines = []
with open(vrtEV_EmissiveDir, 'r') as f:
for line in f:
lines.append(line)
lines.insert(1, '<Metadata domain="GEOLOCATION">\n')
lines.insert(2, ' <MDI key="LINE_OFFSET">1</MDI>\n')
lines.insert(3, ' <MDI key="LINE_STEP">1</MDI>\n')
lines.insert(4, ' <MDI key="PIXEL_OFFSET">1</MDI>\n')
lines.insert(5, ' <MDI key="PIXEL_STEP">1</MDI>\n')
lines.insert(6,
' <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</MDI>\n')
lines.insert(7, ' <MDI key="X_BAND">1</MDI>')
lines.insert(8, ' <MDI key="X_DATASET">HDF5:"{}"://Geolocation/Longitude</MDI>\n'.format(geoFile))
lines.insert(9, ' <MDI key="Y_BAND">1</MDI>\n')
lines.insert(10, ' <MDI key="Y_DATASET">HDF5:"{}"://Geolocation/Latitude</MDI>\n'.format(geoFile))
lines.insert(11, '</Metadata>\n')
with open(vrtEV_EmissiveDir, 'w') as f:
for line in lines:
f.writelines(line)
'''
geoData = gdal.Warp('F:\Py_project\gdalDealMERSI2\Warp_B24_WGS.tif',
'F:\Py_project\gdalDealMERSI2\MERSI2_0403_B24_vrt.vrt',
format='GTiff', geoloc=True, dstSRS="EPSG:4326",
xRes=0.25, yRes=0.25)
'''
# 输出文件名字
os.path.basename(file).split()
tempname1= os.path.splitext(file)[0] + 'EV_Emissive' + '.tif'
tempname2= os.path.splitext(file)[0] + 'EV_RefSB' + '.tif'
EV_Emissivefile = os.path.join(afterGeoPath,os.path.basename(tempname1))
geoData = gdal.Warp(EV_Emissivefile, vrtEV_EmissiveDir,
format='GTiff', geoloc=True, dstSRS=srs,
resampleAlg=gdal.GRIORA_Bilinear, xRes=1000, yRes=1000)
lines = []
with open(vrtEV_RefSBDir, 'r') as f:
for line in f:
lines.append(line)
lines.insert(1, '<Metadata domain="GEOLOCATION">\n')
lines.insert(2, ' <MDI key="LINE_OFFSET">1</MDI>\n')
lines.insert(3, ' <MDI key="LINE_STEP">1</MDI>\n')
lines.insert(4, ' <MDI key="PIXEL_OFFSET">1</MDI>\n')
lines.insert(5, ' <MDI key="PIXEL_STEP">1</MDI>\n')
lines.insert(6,
' <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</MDI>\n')
lines.insert(7, ' <MDI key="X_BAND">1</MDI>')
lines.insert(8, ' <MDI key="X_DATASET">HDF5:"{}"://Geolocation/Longitude</MDI>\n'.format(geoFile))
lines.insert(9, ' <MDI key="Y_BAND">1</MDI>\n')
lines.insert(10, ' <MDI key="Y_DATASET">HDF5:"{}"://Geolocation/Latitude</MDI>\n'.format(geoFile))
lines.insert(11, '</Metadata>\n')
with open(vrtEV_RefSBDir, 'w') as f:
for line in lines:
f.writelines(line)
EV_RefSBfile = os.path.join(afterGeoPath, os.path.basename(tempname2))
geoData = gdal.Warp(EV_RefSBfile, vrtEV_RefSBDir,
format='GTiff', geoloc=True, dstSRS=srs,
resampleAlg=gdal.GRIORA_Bilinear, xRes=1000, yRes=1000)
if geoData == None:
print('deal failure!')
del geoData
print('{} finish Geo\n'.format(file))
if __name__ == '__main__':
afterGeoPath = r'outpath'
fy3file = r'FY3C_VIRRX_GBAL_L1_20211110_0825_1000M_MS.HDF'
fy3GEOfile =r'FY3C_VIRRX_GBAL_L1_20211110_0825_GEOXX_MS.HDF'
geoMERSI2(fy3file, fy3GEOfile, afterGeoPath)
