本人为了做一个气象文件的操作，不知查看了多少文章，问了多少次chatgpt，所以在此本着前人掉坑，后人闪开的想法记录一下。不足之处请留言改正。
# region #主动折叠代码
import os
import xarray as xr
import rioxarray as rio
from datetime import date
import netCDF4 as nc
import rioxarray as rio #!!!这里虽然文件没有主动引用，但是实际上使用了！！！
import matplotlib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import colors
# endregion #折叠结束
#假设nc文件包含经纬度lon,lat和时间time三个刻度
#参数包含 tp降水 sp气压 tm温度
#读取nc文件到内存（xarray.dataset）(文件路径不允许出现中文（毕竟都是是老外发明的工具）)
#单个nc文件
ds = xr.opendataset("file_path.nc")#此处可以使用dask辅助，文章挺多的，不在这里做示范
#多个nc文件按照某个维度合并成一个
file_paths = os.listdir("file_dir")#获取某个目录下所有的文件路径组成列表
'''对象示例
[
"file_path1.nc",
"file_path12.nc",...
]
'''
#此方法会默认使用dask，所以需要安装dask包 否则会报错！
ds = xr.open_mfdataset(file_paths， combine = 'time')
#选中一个参数 tp
da = ds["tp"] #或 da = ds.tp
#如果nc文件只有一个参数 可以直接xr读取
da = xr.open_dataarray('file_path.nc')
#时间选时间刻度的第一个时间
da = da.isel(time=0)
#选择经纬度范围
lon = slice(80, 140)
lat = slice(20, 60)
da = da.sel(lon = lon , lat = lat)
#添加属性
da.attrs["max"] = da.max().values.item()
da.attrs["min"] = da.min().values.item()
#查看指定位置的值
var = da.sel(lon = 100,lat = 30 ,method = 'nearest')
#method='nearest'：选择与目标坐标值最接近的数据。如果存在多个相等距离的坐标值，则选择第一个遇到
#坐标值。
#method='ffill' 或 method='pad'：通过使用前面最近的非缺失值来填充缺失的坐标值。
#method='bfill' 或 method='backfill'：通过使用后面最近的非缺失值来填充缺失的坐标值。
#写入投影方式
da = da.rio.write_crs(4326)
#da = da.rio.write_crs(3856)#一般情况下是4326投影
#4326->3857
#da = da.rio.write_crs(4326).rio.reproject(
#"EPSG:3857"
#).rio.write_crs(3857)
#注意:3857的裁剪需要坐标转换器
#crs_3857 = pyproj.CRS.from_string("+init=epsg:3857")
#lonlat_to_3857 = pyproj.Transformer.from_crs(
"EPSG:4326", crs_3857, always_xy=True
)
#x_min, y_min = lonlat_to_3857.transform(80, 140)
#x_max, y_max = lonlat_to_3857.transform(20, 60)
#da = da.rio.clip_box(x_min, y_min, x_max, y_max)
#设置分辨率0.1(分辨率越小，画质越高)
#注意3857的分辨率在4326的时候设置，转化成3857后执行此代码报错
da = da.rio.reproject(
"EPSG:4326",
resolution=0.1,
)
#插入空值
da = da.rio.interpolate_na(method='linear')
#插值方法可以选择多种：
#线性插值、最近邻插值或三次样条插值
#生成tif,tiff文件
da.rio.to_raster("tif_name.tif")
下面画图
#包和数据是上面的
#调制色标(颜色支持rgba(1,1,1,1)与十六进制颜色百分比透明度)
color_list = ["red", "orange", "yellow", "green","Cyan", "blue", "purple","black"]
cmap = colors.ListedColormap(color_list )
#对应的数据
bounds = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1]
#适应气象数据值
var = da.attr["max"] - da.attr["min"]
bounds = [x * var + da.attr["min"] for x in bounds]
norm = colors.BoundaryNorm(bounds, cmap.N)
color = {
"cmap": cmap,
"norm": norm,
}
# 创建画布和坐标轴
fig, ax = plt.subplots(
figsize=(10, 8),#图片大小
subplot_kw={"projection": ccrs.PlateCarree()}#地图投影
)
da.plot(
**color,#使用自定义颜色 也可以使用内置 cmap = "jet"等等
ax=ax,
transform=ccrs.PlateCarree(),
cbar_kwargs={"label": None, "location": "bottom"},
)
# 添加河流和国界
ax.add_feature(cfeature.RIVERS)#河流
ax.add_feature(cfeature.BORDERS)#国界
ax.add_feature(cfeature.COASTLINE)#海岸线
ax.add_feature(cfeature.LAND)#陆地
ax.add_feature(cfeature.LAKES)#湖泊
provinces = cfeature.NaturalEarthFeature(
category="cultural",
name="admin_1_states_provinces_lines",
scale="50m",
facecolor="none",
)
ax.add_feature(provinces, edgecolor="black", linewidth=0.5)#省界线
# 设置x,y轴标签
ax.set_xticks(da.lon.values.tolist(), "lon")
ax.set_yticks(da.lat.values.tolist(), "lat")
# 坐标轴名称
ax.set_title("tp", size=20)#默认不支持中文，可以设置，有些麻烦，不想在这写
ax.set_xlabel("lon")
ax.set_ylabel("lat")
#展示
plt.show()
#保存成png图片
fig.savefig("img.png", dpi=100)#dpi设置图片分辨率:分辨率越大，图片越精细，文件越大!
就写这些吧，最简单的nc处理和图片展示
文章链接
您阅读本篇文章共花了：
发表评论