Контурная схема в Python для вывода нескольких netcdf - PullRequest
Новая
       14

Контурная схема в Python для вывода нескольких netcdf

0 голосов
/ 13 марта 2019

В настоящее время я строю карту ветра, основанную на месяце, используя wrf-python.Я столкнулся с трудностями в plt.contourf().Мои lons (x) и lats (y) имеют 2D (south_north, west_east), в то время как мои wspd_500 (z) имеют 3D (time, south_north, west_east), так как я вызвал 2 файла wrfoutput.Вот почему, когда я запускаю скрипт, ошибка говорит о том, что z должно быть 2D.Как мне решить это? 

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature

from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy,
                 cartopy_xlim, cartopy_ylim, ALL_TIMES)



# Open the NetCDF file
ncfile = [Dataset("wrfout_d01_2016-08-01_01:00:00"),
          Dataset("wrfout_d01_2016-08-01_01:00:00")]

# Extract the pressure, geopotential height, and wind variables
p = getvar(ncfile, "pressure")
z = getvar(ncfile, "z", units="dm")
#ua = getvar(ncfile, "ua", timeidx=ALL_TIMES, method="join", units="m s-1")
#va = getvar(ncfile, "va", timeidx=ALL_TIMES, method="join", units="m s-1")
wspd = getvar(ncfile, "wspd_wdir", timeidx=ALL_TIMES, method="cat", units="m s-1")[0,:]

# Interpolate geopotential height, u, and v winds to 500 hPa
ht = interplevel(z, p, 500)  
#u_500 = interplevel(ua, z, 500)  
#v_500 = interplevel(va, z, 500)  
wspd_500 = interplevel(wspd, z, 500)  

# Get the lat/lon coordinates
lats, lons = latlon_coords(ht)

print(wspd_500.dims)

# Get the map projection information
cart_proj = get_cartopy(ht)

# Create the figure
fig = plt.figure(figsize=(12,9))
ax = plt.axes(projection=cart_proj)

# Download and add the states and coastlines
states = NaturalEarthFeature(category="cultural", scale="50m",
                             facecolor="none",
                             name="admin_1_states_provinces_shp")
ax.add_feature(states, linewidth=1.0, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)

# Add the wind speed contours
levels = np.arange(1., 20., 2.)
wspd_contours = plt.contourf(to_np(lons), to_np(lats), to_np(wspd_500),
                             levels=levels,
                             cmap=get_cmap("rainbow"),
                             transform=crs.Mercator())
plt.colorbar(wspd_contours, ax=ax, orientation="horizontal", pad=.05)

# Set the map bounds
ax.set_xlim(cartopy_xlim(ht))
ax.set_ylim(cartopy_ylim(ht))

ax.gridlines()

plt.title("Wind Speed (m s-1), Barbs (kt)")

plt.show()
...