格点数据可视化(美国站点的日降雨数据)

获取美国站点的日降雨量的格点数据,并且可视化
在这里插入图片描述

导入模块

from datetime import datetime, timedelta
from urllib.request import urlopen

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from metpy.units import masked_array, units
from netCDF4 import Dataset

读取数据

nc = Dataset('20200309_conus.nc')
prcpvar = nc.variables['observation']
data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('mm')
x = nc.variables['x'][:]
y = nc.variables['y'][:]
proj_var = nc.variables[prcpvar.grid_mapping]

设置投影

globe = ccrs.Globe(semimajor_axis=proj_var.earth_radius)
proj = ccrs.Stereographic(central_latitude=90.0,
                          central_longitude=proj_var.straight_vertical_longitude_from_pole,
                          true_scale_latitude=proj_var.standard_parallel, globe=globe)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)

# 绘制海岸线、国界线、州界线
ax.coastlines()
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.STATES)

# 设置降雨量等级间隔
clevs = [0, 1, 2.5, 5, 7.5, 10, 15, 20, 30, 40,
         50, 70, 100, 150, 200, 250, 300, 400, 500, 600, 750]
# In future MetPy
# norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)
# 单独设置cmap
cmap_data = [(1.0, 1.0, 1.0),
             (0.3137255012989044, 0.8156862854957581, 0.8156862854957581),
             (0.0, 1.0, 1.0),
             (0.0, 0.8784313797950745, 0.501960813999176),
             (0.0, 0.7529411911964417, 0.0),
             (0.501960813999176, 0.8784313797950745, 0.0),
             (1.0, 1.0, 0.0),
             (1.0, 0.6274510025978088, 0.0),
             (1.0, 0.0, 0.0),
             (1.0, 0.125490203499794, 0.501960813999176),
             (0.9411764740943909, 0.250980406999588, 1.0),
             (0.501960813999176, 0.125490203499794, 1.0),
             (0.250980406999588, 0.250980406999588, 1.0),
             (0.125490203499794, 0.125490203499794, 0.501960813999176),
             (0.125490203499794, 0.125490203499794, 0.125490203499794),
             (0.501960813999176, 0.501960813999176, 0.501960813999176),
             (0.8784313797950745, 0.8784313797950745, 0.8784313797950745),
             (0.9333333373069763, 0.8313725590705872, 0.7372549176216125),
             (0.8549019694328308, 0.6509804129600525, 0.47058823704719543),
             (0.6274510025978088, 0.42352941632270813, 0.23529411852359772),
             (0.4000000059604645, 0.20000000298023224, 0.0)]
            
cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)

cs = ax.contourf(x, y, data, clevs, cmap=cmap, norm=norm)

# 添加colorbar
cbar = plt.colorbar(cs, orientation='horizontal')
cbar.set_label(data.units)
# 设置标题
ax.set_title(prcpvar.long_name + ' for period ending ' + nc.creation_time)
plt.show()

数据怎样获取

dt = datetime.utcnow() - timedelta(days=1)  # 获取过去1天的时间
url = ('http://water.weather.gov/precip/downloads/{dt:%Y/%m/%d}/nws_precip_1day_'
       '{dt:%Y%m%d}_conus.nc'.format(dt=dt))
data = urlopen(url).read()
nc = Dataset('data', memory=data)

显示数据

import xarray as xr
from xarray.backends import NetCDF4DataStore
data = xr.open_dataset(NetCDF4DataStore(nc))
data

保存为nc数据文章来源地址https://www.uudwc.com/A/dbxOo/

data.to_netcdf('{dt:%Y%m%d}_conus.nc'.format(dt=dt),'w')

原文地址:https://blog.csdn.net/weixin_45492560/article/details/133441978

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请联系站长进行投诉反馈,一经查实,立即删除!

h
上一篇 2023年10月24日 08:04
spring6概述
下一篇 2023年10月24日 09:04