VNP46A3灯光遥感数据全球拼接并重采样

感谢Deepseek帮我写代码，本人在此过程中仅对其进行调试和部分修改：
灯光遥感2024年1月全球拼接结果
代码如下：

import os
import glob
import h5py
import numpy as np
from osgeo import gdal, osr
import rasterio
from rasterio.merge import merge
from rasterio.warp import reproject, Resampling, calculate_default_transform
import tempfiledef resample_raster(input_path, output_path, target_resolution=0.5):"""重采样栅格数据并正确处理无效值"""with rasterio.open(input_path) as src:# 读取原始nodata值（如果没有则使用65535）src_nodata = 65535# 计算目标transform和尺寸dst_crs = src.crs  # 保持原坐标系transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds, resolution=(target_resolution, target_resolution))# 设置输出metadatadst_meta = src.meta.copy()dst_meta.update({'transform': transform,'width': width,'height': height,'nodata': src_nodata  # 保持一致的nodata值})# 执行重采样（自动排除无效值）with rasterio.open(output_path, 'w', **dst_meta) as dst:for i in range(1, src.count + 1):# 方法1：使用rasterio的自动nodata处理reproject(source=rasterio.band(src, i),destination=rasterio.band(dst, i),src_transform=src.transform,src_crs=src.crs,dst_transform=transform,dst_crs=dst_crs,resampling=Resampling.average,src_nodata=src_nodata,  # 关键：指定源无效值dst_nodata=src_nodata   # 保持输出无效值一致)def process_viirs_to_global_tiff(input_folder, output_path, target_resolution=0.5):"""处理VIIRS数据并直接重采样到指定分辨率"""# 获取所有h5文件h5_files = glob.glob(os.path.join(input_folder, '*.h5'))# 临时存储所有分块GeoTIFFtemp_files = []# 处理每个h5文件for h5_file in h5_files:try:with h5py.File(h5_file, 'r') as f:# 读取数据data = f['HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/AllAngle_Composite_Snow_Free'][:]lon = f['HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/lon'][:]lat = f['HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/lat'][:]# data[data == 65535] = np.nan    # 将无效值设为NaN# 创建临时GeoTIFFtemp_file = tempfile.NamedTemporaryFile(suffix='.tif', delete=False).nametemp_files.append(temp_file)# 获取数据维度rows, cols = data.shape# 创建输出GeoTIFFdriver = gdal.GetDriverByName('GTiff')out_ds = driver.Create(temp_file, cols, rows, 1, gdal.GDT_Float32)# 设置地理参考信息（使用正弦投影）srs = osr.SpatialReference()srs.ImportFromEPSG(4326)# srs.ImportFromProj4("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")out_ds.SetProjection(srs.ExportToWkt())# 计算地理变换min_lon, max_lon = np.min(lon), np.max(lon)min_lat, max_lat = np.min(lat), np.max(lat)pixel_width = (max_lon - min_lon) / colspixel_height = (max_lat - min_lat) / rowsgeo_transform = (min_lon, pixel_width, 0, max_lat, 0, -pixel_height)out_ds.SetGeoTransform(geo_transform)# 写入数据out_band = out_ds.GetRasterBand(1)out_band.WriteArray(data)out_band.FlushCache()# out_band.SetNoDataValue(-9999)out_ds = None  # 关闭文件except Exception as e:print(f"处理文件 {h5_file} 时出错: {str(e)}")continue# 合并所有临时GeoTIFFif temp_files:try:# 创建临时合并文件merged_temp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False).name# 打开所有临时文件src_files_to_mosaic = []for file in temp_files:src = rasterio.open(file)src_files_to_mosaic.append(src)# 合并mosaic, out_trans = merge(src_files_to_mosaic)# 写入临时合并文件out_meta = src_files_to_mosaic[0].meta.copy()out_meta.update({"driver": "GTiff","height": mosaic.shape[1],"width": mosaic.shape[2],"transform": out_trans})with rasterio.open(merged_temp, "w", **out_meta) as dest:dest.write(mosaic)# 对合并后的文件进行重采样resample_raster(merged_temp, output_path, target_resolution)print(f"成功创建并重采样全球TIFF文件: {output_path}")except Exception as e:print(f"合并或重采样文件时出错: {str(e)}")finally:# 关闭所有文件并清理临时文件for src in src_files_to_mosaic:src.close()for file in temp_files + [merged_temp]:try:os.remove(file)except:passelse:print("没有有效的HDF5文件处理完成")# 使用示例
input_folder = r"E:\DBN_data\VNP46A3\2024\001"
day = int(input_folder[-3:])
output_tiff = os.path.join(input_folder, 'DNB_merge_'+str(input_folder[-8:-4])+str(day).zfill(2)+'_05deg.tif')
process_viirs_to_global_tiff(input_folder, output_tiff, target_resolution=0.5)

其中用的较多的，rasterio库，这个库我自己确实不知道这么多用处，感谢AI，让我学到了新东西。