问题
马来西亚的合作伙伴咨询说CASA提供的Tansat XCO2数据无法正常读取和转换,我尝试写了一段代码读取国家对地观测科学数据中心网站上发布的V2版本的逐日数据,验证CASA的数据没有问题。
思路
读取Tansat XCO2的nc文件 获取 经纬度和XCO2浓度,然后转换为shp点数据,输出。
有的nc文件直接用arcgis的Make NetCDF Raster Layer功能就可以直接转换输出,有的自变量不是经纬度,经纬度也存在numpy的数组中,需要提取来。
步骤
运行环境:Python 3.11.5原始的jupyter notebook文件点击此处下载。
1、读取文件
import netCDF4 as nc
import pandas as pd
import geopandas as gpd
#读取 .nc文件
#TanSat_L2_XCO2_IAPCAS_V2 (Daily data)
file_path = r'M:\WT\Tansat\TanSat_L2_XCO2_IAPCAS_V2\daily\TanSat_ACGS_SCI_ND_L2_XCO2_lite_20170301.nc'
dataset = nc.Dataset(file_path)
2、查看经纬度和XCO2变量的存储列名
#获得变量列
variables = dataset.variables
variables
3、 创建DataFrame,获取经纬度、XCO2浓度和时间
df_xco2=pd.DataFrame()
df_xco2['xco2']=dataset.variables['xco2'][:]
df_xco2['longitude']=dataset.variables['longitude'][:]
df_xco2['latitude']=dataset.variables['latitude'][:]
df_xco2['datetime']=dataset.variables['sounding_id'][:]
df_xco2
4、 将DataFrame转换为GeoDataFrame,而后转为shp点数据输出
# 创建GeoDataFrame,导出shp数据
gdf = gpd.GeoDataFrame(
df_xco2,
geometry=gpd.points_from_xy(df_xco2.longitude, df_xco2.latitude)
)
# 设置坐标参考系统为WGS 84
gdf.set_crs(epsg=4326, inplace=True)
# 保存为Shapefile
gdf.to_file("TanSat_ACGS_SCI_ND_L2_XCO2_lite_20170301.shp")
5、在QGIS中加载shp数据
引用文献:
[0] ArcGIS:读取nc格式文件并导出为tif格式文件,降雨或温度NC等数据