MODIS以及VIIRS的一级产品以及部分的大气产品都是通过Swath的方式进行存储,即每个像元对应的面积大小是不一致的。一般这种数据存储都是每个像元都会对应一个经纬度坐标,不能通过ArcGIS进行处理。 如何把Swath数据转化成每个像元面积相等的,即Grid数据就是一个比较基础数据预处理过程。常用的方法就是用ENVI、IDL之类的进行处理。今天我分享使用Python中的pyresample库对Swath数据进行几何校正。
样例数据使用的是VIIRS的LST产品,产品分辨率为750m,产品名为:VNP21,数据为:VNP21.A2022365.0624.001.2023020015248.nc。
未进行几何校正的数据
可以看出直接把这个数据放到ArcGIS里面,已经偏的不行了。这景数据的实际覆盖范围其实是我国的西南部以及部分东南亚区域。
下面是我使用Python对数据进行几何校正的代码,大家可以参考一下,如果有什么问题,也希望能够及时指正。
# -*- coding: utf-8 -*-
from netCDF4 import Dataset
from pyproj import CRS
import numpy as np
from pyresample import geometry as geom
from pyresample import kd_tree as kdt
from osgeo import gdal, osr
def h5_tif(data_path):
f = Dataset(data_path, 'r')
Data_Fields_variables=['Emis_14','Emis_15','Emis_16','Emis_ASTER','LST']
for Data_Fields_variable in Data_Fields_variables:
Data_Fields_data=f['VIIRS_Swath_LSTE']['Data Fields'][Data_Fields_variable][:]
Data_Fields_scale=f['VIIRS_Swath_LSTE']['Data Fields'][Data_Fields_variable].scale_factor
Data_Fields_offset=f['VIIRS_Swath_LSTE']['Data Fields'][Data_Fields_variable].add_offset
Data_Fields_FillValue=f['VIIRS_Swath_LSTE']['Data Fields'][Data_Fields_variable]._FillValue
Data_Fields_data=np.where(Data_Fields_data!=Data_Fields_FillValue,
Data_Fields_data*Data_Fields_scale+Data_Fields_offset,
Data_Fields_FillValue)
lon=f['VIIRS_Swath_LSTE']['Geolocation Fields']['longitude'][:]
lat=f['VIIRS_Swath_LSTE']['Geolocation Fields']['latitude'][:]
lon[lon==f['VIIRS_Swath_LSTE']['Geolocation Fields']['longitude']._FillValue]=np.nan
lat[lat==f['VIIRS_Swath_LSTE']['Geolocation Fields']['latitude']._FillValue]=np.nan
epsg, proj, pName = '4326', 'longlat', 'Geographic' # Set output projection to Geographic CRS
llLon, llLat, urLon, urLat = np.nanmin(lon), np.nanmin(lat), np.nanmax(lon), np.nanmax(lat)
areaExtent = (llLon, llLat, urLon, urLat)
ps = 0.02 # Hard-coded estimate of pixel size at 0.02 degrees
cols = int(round((areaExtent[2] - areaExtent[0]) / ps)) # Calculate the output cols
rows = int(round((areaExtent[3] - areaExtent[1]) / ps)) # Calculate the output rows
areaDef = geom.AreaDefinition(epsg, pName, proj, CRS.from_epsg(4326), cols, rows, areaExtent)
swath_def = geom.SwathDefinition(lons=lon, lats=lat)
result = kdt.resample_nearest(swath_def, Data_Fields_data,areaDef, radius_of_influence=5000)
gt = [areaDef.area_extent[0], ps, 0, areaDef.area_extent[3], 0, -ps]
# Get driver, specify dimensions, define and set output geotransform
height, width = result.shape # Define geotiff dimensions
driver = gdal.GetDriverByName('GTiff')
d = driver.Create(data_path.replace('.nc',Data_Fields_variable+'.tif'), width, height, 1, gdal.GDT_Float32)
d.SetGeoTransform(gt)
srs = osr.SpatialReference()
srs.ImportFromEPSG(int(epsg))
d.SetProjection(srs.ExportToWkt())
band = d.GetRasterBand(1)
band.WriteArray(result)
band.SetNoDataValue(int(Data_Fields_FillValue))
band.FlushCache()
if __name__ == "__main__":
h5_tif(r'F:\code\swath2grid\VNP21.A2022365.0624.001.2023020015248.nc')
几何校正的结果
这就是几何校正后的结果,可以看出来这个结果还是蛮不错的,空间位置与shp还是非常一致的。