Python处理NetCDF格式数据（以PM2.5数据为例）并导出为TIFF数据(附源代码)

最新推荐文章于 2024-01-26 12:15:13 发布

我要抓住你啦

最新推荐文章于 2024-01-26 12:15:13 发布

阅读量1.3k

点赞数 2

分类专栏： Python处理PM2.5 文章标签： python

本文链接：https://blog.csdn.net/weixin_38605438/article/details/110500393

版权

Python处理PM2.5 专栏收录该内容

1 篇文章 1 订阅

订阅专栏

Python处理NetCDF格式数据（以PM2.5数据为例）并导出为TIFF数据(附源代码)

Python处理NetCDF格式数据

以PM2.5数据为例

PM2.5数据出处：Atmospheric Composition Analysis Group

并导出为TIFF

-------------------------------------------------------

所用代码如下：

存在一个不太理解问题，还请各位指教：
读取原始nc数据后，需将提取的PM2.5逆时针旋转90°才可以

#!usr/bin/env python
# -*- coding: utf-8 -*-

import os
import gdal,osr
import netCDF4 as nc
import numpy as np

def	extract_nc(nc_file):
	dataset = nc.Dataset(nc_file)
	print(dataset.variables.keys())
	
	# 读取经纬度
	var_lat = dataset.variables["LAT"][:]
	var_lon = dataset.variables["LON"][:]

	PM25_1 = dataset.variables["PM25"]
	PM25 = np.rot90(PM25_1,1) # 此处我没有找到比较好的解决方法，若不进行此处的逆时针旋转，数据集输出将会有误，且会被旋转90°

	# 影像的左上角和右下角坐标
	LonMin, LatMax, LonMax, LatMin = [var_lon.min(), var_lat.max(), var_lon.max(), var_lat.min()]
	# 分辨率计算，可以计算得到
	# 也可以自己手动输入，我所用数据集的分辨率是0.01°×0.01°
	N_Lat = len(var_lat)
	N_Lon = len(var_lon)
	Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
	Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
 
	# 设置影像的显示范围
	# Lat_Res一定要是-的，此处取负号原因不是很明白
	geo = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)

	#构造projection
	src_srs = osr.SpatialReference()
	src_srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
	src_srs_wkt = src_srs.ExportToWkt()
	return PM25,geo,src_srs_wkt,N_Lon,N_Lat

def write_tiff(tiff_file,geo,proj,N_Lon,N_Lat,data_array):
	driver = gdal.GetDriverByName("Gtiff")

	outdataset = driver.Create(tiff_file,N_Lon,N_Lat,1,gdal.GDT_Float32)

	outdataset.SetGeoTransform(geo)
	outdataset.SetProjection(proj)

	band = outdataset.GetRasterBand(1)
	band.WriteArray(data_array)

def nc_to_tif(nc_file,output_dir):

	nc_file_name = os.path.basename(nc_file)
	output_file_name = os.path.join(output_dir,nc_file_name)

	pm25_file = output_file_name.replace(".nc", "_2000.tif")

	pm25, geo, proj,N_Lon,N_Lat=extract_nc(nc_file)

	write_tiff(pm25_file,geo,proj,N_Lon,N_Lat,pm25)

if __name__ == '__main__':
	# PM2.5数据
	PM_file = r"E:\AYiDaiYiLu\PM\2000PM.nc"
	# 解析结果存储路径
	tiff_dir = r"E:\AYiDaiYiLu"
	nc_to_tif(PM_file,tiff_dir)