当前位置:   article > 正文

python气象(地理)nc数据处理 使用xarray matplotlib rioxarray库完成nc文件转tif(tiff)文件,绘制地理图片_xarray保存nc文件

xarray保存nc文件

本人为了做一个气象文件的操作,不知查看了多少文章,问了多少次chatgpt,所以在此本着前人掉坑,后人闪开的想法记录一下。不足之处请留言改正。

  1. # region #主动折叠代码
  2. import os
  3. import xarray as xr
  4. import rioxarray as rio
  5. from datetime import date
  6. import netCDF4 as nc
  7. import rioxarray as rio #!!!这里虽然文件没有主动引用,但是实际上使用了!!!
  8. import matplotlib
  9. import matplotlib.pyplot as plt
  10. import cartopy.crs as ccrs
  11. import cartopy.feature as cfeature
  12. from matplotlib import colors
  13. # endregion #折叠结束
  14. #假设nc文件包含经纬度lon,lat和时间time三个刻度
  15. #参数包含 tp降水 sp气压 tm温度
  16. #读取nc文件到内存(xarray.dataset)(文件路径不允许出现中文(毕竟都是是老外发明的工具))
  17. #单个nc文件
  18. ds = xr.opendataset("file_path.nc")#此处可以使用dask辅助,文章挺多的,不在这里做示范
  19. #多个nc文件按照某个维度合并成一个
  20. file_paths = os.listdir("file_dir")#获取某个目录下所有的文件路径组成列表
  21. '''对象示例
  22. [
  23. "file_path1.nc",
  24. "file_path12.nc",...
  25. ]
  26. '''
  27. #此方法会默认使用dask,所以需要安装dask包 否则会报错!
  28. ds = xr.open_mfdataset(file_paths, combine = 'time')
  29. #选中一个参数 tp
  30. da = ds["tp"] #或 da = ds.tp
  31. #如果nc文件只有一个参数 可以直接xr读取
  32. da = xr.open_dataarray('file_path.nc')
  33. #时间选时间刻度的第一个时间
  34. da = da.isel(time=0)
  35. #选择经纬度范围
  36. lon = slice(80, 140)
  37. lat = slice(20, 60)
  38. da = da.sel(lon = lon , lat = lat)
  39. #添加属性
  40. da.attrs["max"] = da.max().values.item()
  41. da.attrs["min"] = da.min().values.item()
  42. #查看指定位置的值
  43. var = da.sel(lon = 100,lat = 30 ,method = 'nearest')
  44. #method='nearest':选择与目标坐标值最接近的数据。如果存在多个相等距离的坐标值,则选择第一个遇到
  45. #坐标值。
  46. #method='ffill' 或 method='pad':通过使用前面最近的非缺失值来填充缺失的坐标值。
  47. #method='bfill' 或 method='backfill':通过使用后面最近的非缺失值来填充缺失的坐标值。
  48. #写入投影方式
  49. da = da.rio.write_crs(4326)
  50. #da = da.rio.write_crs(3856)#一般情况下是4326投影
  51. #4326->3857
  52. #da = da.rio.write_crs(4326).rio.reproject(
  53. #"EPSG:3857"
  54. #).rio.write_crs(3857)
  55. #注意:3857的裁剪需要坐标转换器
  56. #crs_3857 = pyproj.CRS.from_string("+init=epsg:3857")
  57. #lonlat_to_3857 = pyproj.Transformer.from_crs(
  58. "EPSG:4326", crs_3857, always_xy=True
  59. )
  60. #x_min, y_min = lonlat_to_3857.transform(80, 140)
  61. #x_max, y_max = lonlat_to_3857.transform(20, 60)
  62. #da = da.rio.clip_box(x_min, y_min, x_max, y_max)
  63. #设置分辨率0.1(分辨率越小,画质越高)
  64. #注意3857的分辨率在4326的时候设置,转化成3857后执行此代码报错
  65. da = da.rio.reproject(
  66. "EPSG:4326",
  67. resolution=0.1,
  68. )
  69. #插入空值
  70. da = da.rio.interpolate_na(method='linear')
  71. #插值方法可以选择多种:
  72. #线性插值、最近邻插值或三次样条插值
  73. #生成tif,tiff文件
  74. da.rio.to_raster("tif_name.tif")

下面画图

  1. #包和数据是上面的
  2. #调制色标(颜色支持rgba(1,1,1,1)与十六进制颜色百分比透明度)
  3. color_list = ["red", "orange", "yellow", "green","Cyan", "blue", "purple","black"]
  4. cmap = colors.ListedColormap(color_list )
  5. #对应的数据
  6. bounds = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1]
  7. #适应气象数据值
  8. var = da.attr["max"] - da.attr["min"]
  9. bounds = [x * var + da.attr["min"] for x in bounds]
  10. norm = colors.BoundaryNorm(bounds, cmap.N)
  11. color = {
  12. "cmap": cmap,
  13. "norm": norm,
  14. }
  15. # 创建画布和坐标轴
  16. fig, ax = plt.subplots(
  17. figsize=(10, 8),#图片大小
  18. subplot_kw={"projection": ccrs.PlateCarree()}#地图投影
  19. )
  20. da.plot(
  21. **color,#使用自定义颜色 也可以使用内置 cmap = "jet"等等
  22. ax=ax,
  23. transform=ccrs.PlateCarree(),
  24. cbar_kwargs={"label": None, "location": "bottom"},
  25. )
  26. # 添加河流和国界
  27. ax.add_feature(cfeature.RIVERS)#河流
  28. ax.add_feature(cfeature.BORDERS)#国界
  29. ax.add_feature(cfeature.COASTLINE)#海岸线
  30. ax.add_feature(cfeature.LAND)#陆地
  31. ax.add_feature(cfeature.LAKES)#湖泊
  32. provinces = cfeature.NaturalEarthFeature(
  33. category="cultural",
  34. name="admin_1_states_provinces_lines",
  35. scale="50m",
  36. facecolor="none",
  37. )
  38. ax.add_feature(provinces, edgecolor="black", linewidth=0.5)#省界线
  39. # 设置x,y轴标签
  40. ax.set_xticks(da.lon.values.tolist(), "lon")
  41. ax.set_yticks(da.lat.values.tolist(), "lat")
  42. # 坐标轴名称
  43. ax.set_title("tp", size=20)#默认不支持中文,可以设置,有些麻烦,不想在这写
  44. ax.set_xlabel("lon")
  45. ax.set_ylabel("lat")
  46. #展示
  47. plt.show()
  48. #保存成png图片
  49. fig.savefig("img.png", dpi=100)#dpi设置图片分辨率:分辨率越大,图片越精细,文件越大!

就写这些吧,最简单的nc处理和图片展示

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/羊村懒王/article/detail/329860
推荐阅读
相关标签
  

闽ICP备14008679号