0

When I run this code 

import Scientific.IO.NetCDF as S
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import xarray as xr
import metpy
import numpy as N

from metpy.plots import ContourPlot, ImagePlot, MapPanel, PanelContainer
# Any import of metpy will activate the accessors
import metpy.calc as mpcalc
#from metpy.testing import get_test_data
from metpy.units import units
# Open the netCDF file as a xarray Datase


#
datadir='C:/Users/stratus/AppData/Local/lxss/home/stratus/PROJECT/NEWPROJECT/FEB012017/nam_218_20170131_1200_000.nc'
data = xr.open_dataset(datadir,decode_cf=True)
# To parse the full dataset, we can call parse_cf without an argument, and assign the returned
# Dataset.
data = data.metpy.parse_cf()
tempatt=data['TMP_P0_L100_GLC0'].attrs

# If we instead want just a single variable, we can pass that variable name to parse_cf and
# it will return just that data variable as a DataArray.
data_var = data.metpy.parse_cf('TMP_P0_L100_GLC0')

# To rename variables, supply a dictionary between old and new names to the rename method
data.rename({
    'TMP_P0_L100_GLC0': 'temperature',
}, inplace=True)

data['temperature'].metpy.convert_units('degC')

# Get multiple coordinates (for example, in just the x and y direction)
x, y = data['temperature'].metpy.coordinates('x', 'y')

# If we want to get just a single coordinate from the coordinates method, we have to use
# tuple unpacking because the coordinates method returns a generator
vertical, = data['temperature'].metpy.coordinates('vertical')
data_crs = data['temperature'].metpy.cartopy_crs


# Or, we can just get a coordinate from the property
#time = data['temperature'].metpy.time

# To verify, we can inspect all their names
#print([coord.name for coord in (x, y, vertical, time)])

#
#heights = data['height'].metpy.loc[{'time': time[0], 'vertical': 850. * units.hPa}]
#lat, lon = xr.broadcast(y, x)
#f = mpcalc.coriolis_parameter(lat)
#dx, dy = mpcalc.grid_deltas_from_dataarray(heights)
#u_geo, v_geo = mpcalc.geostrophic_wind(heights, f, dx, dy)
#print(u_geo)
#print(v_geo)
fig=plt.figure(1)
# A very simple example example of a plot of 500 hPa heights

data_crs = data['temperature'].metpy.cartopy_crs

ax = plt.axes(projection=ccrs.LambertConformal())
data['temperature'].metpy.loc[{'vertical': 850. * units.hPa}].plot(ax=ax, transform=data_crs)

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
plt.show()
#ax.set_extent([-120,-80,20,50])
plt.title("850 mb Temperature")
#plt.suptitle("Metpy Test")
plt.show()

I had to edit the code as per some of the answers but I am getting a mostly blank map now. 850 T Map fail I am mainly trying to have the temperatures at 850 mb overlap the US so I could show it to a friend to practice for a project I am helping him with. The filling of the parentheses for the data helped a bit which is why I edited it.

stratusshow
  • 25
  • 9
  • 1
    Any example here (https://unidata.github.io/MetPy/latest/examples/index.html#plotting) that would help? Also, does it change something if you set `transform=ccrs.PlateCarree()` in your call to `plot`? And I think it is best practice to first define Figure and Axes objects before calling `plot`. – Patol75 Oct 11 '19 at 02:42
  • The transform did nothing to the plot. I mentioned that setting the axes came up with a blank png – stratusshow Oct 11 '19 at 05:26
  • 2
    Would it be possible for you to supply data and metpy commands in order to be able to reproduce your map? Then it would be easier to figure out what is going wrong. – Patol75 Oct 11 '19 at 05:53
  • 1
    A full example script would be great to help! I know that where you have plot called won't work since you haven't passed transformation information to it. You would need something like `ax=ax, transform=ccrs.PlateCaree()` inside the `plot()` command. The Xarray tutorial demonstrates some of this [here](https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#plotting). – zbruick Oct 11 '19 at 14:27

1 Answers1

0

As pointed out in the comments it is difficult to answer without a reproducible example. However, the following may solve your issue:

data_crs = data['temperature'].metpy.cartopy_crs

ax = plt.axes(projection=ccrs.LambertConformal())
data['temperature'].metpy.loc[{'vertical': 1000. * units.hPa}].plot(ax=ax, transform=data_crs)

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
plt.show()
baccandr
  • 1,090
  • 8
  • 15
  • I made the modifications and edited the code as such but it still isn't giving me the temperatures over the US like I am trying to intend it for. – stratusshow Oct 11 '19 at 20:51
  • Can you provide the dataset so that we can try to plot the map? By the way you should clean up your code and try to make a more clear minimal, reproducible example: https://stackoverflow.com/help/minimal-reproducible-example – baccandr Oct 14 '19 at 10:12