提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  后端

当前栏目

Swath数据几何校正（Python）

Python数据 几何 校正
2023-06-13 09:15:38 时间

Swath数据几何校正

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还是非常一致的。

相关文章