赞
踩
本人为了做一个气象文件的操作,不知查看了多少文章,问了多少次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处理和图片展示
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。