一种特殊的NC文件转TIF——知道每个像元的坐标值怎么转为TIF

 

 

 

 

import os
import numpy as np
import netCDF4 as nc
from osgeo import gdal

ds = nc.Dataset('./0228SO2.nc')
data = ds['PRODUCT']['sulfurdioxide_total_vertical_column'][0, ...].filled(np.NAN)

lon = ds['PRODUCT']['longitude'][0, ...].data
lat = ds['PRODUCT']['latitude'][0, ...].data

with open('./0228SO2.csv', 'w') as fp:
    fp.write('x,y,z\n')
    for x, y, z in zip(lon.flatten(), lat.flatten(), data.flatten()):
        if not np.isnan(z):
            fp.write(f'{x},{y},{z}\n')

vrt_path = './0228SO2.vrt'
with open(vrt_path, 'w') as fn_vrt:
    layer_name = '0228SO2'
    fn_vrt.write('<OGRVRTDataSource>\n')
    fn_vrt.write('\t<OGRVRTLayer name="%s">\n' % layer_name)
    # 下面这个路径最好使用绝对路径
    fn_vrt.write('\t\t<SrcDataSource>%s</SrcDataSource>\n' % './0228SO2.csv')
    fn_vrt.write('\t\t<GeometryType>wkbPoint</GeometryType>\n')
    fn_vrt.write('\t\t<GeometryField encoding="PointFromColumns" x="x" y="y" z="z"/>\n')
    fn_vrt.write('\t</OGRVRTLayer>\n')
    fn_vrt.write('</OGRVRTDataSource>\n')

# 命令中的各个命令最好都用绝对路径
# 自行设置输出的分辨率或者输出的行列数，以及输入的坐标系
tif_path = './0228SO2.tif'
order = f'gdal_grid -a invdist:power=2.0:smoothing=1.0:nodata=-10000 -txe {lon.min()} {lon.max()} -tye {lat.min()} {lat.max()} -outsize 594 416 -a_srs EPSG:4326 -l {layer_name} {vrt_path} {tif_path} --config GDAL_NUM_THREADS ALL_CPUS'
result = os.system(order)
print(result)

 

 

链接: https://pan.baidu.com/s/1bcBalOLX9HUUtjQ66C-vwQ?pwd=htm7

提取码: htm7