Data with projected coordinates using Xarray and Cartopy

In [1]:

    

%matplotlib inline
import xarray
import matplotlib.pyplot as plt
import numpy as np


    

In [2]:

    

url = 'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/HRRR/CONUS_2p5km/Best'


    

In [3]:

    

nc = xarray.open_dataset(url)


    

In [4]:

    

var='Temperature_height_above_ground'
ncvar = nc[var]
ncvar


    

    
Out[4]:





<xarray.DataArray 'Temperature_height_above_ground' (time: 84, height_above_ground1: 1, y: 1377, x: 2145)>
[248107860 values with dtype=float32]
Coordinates:
  * x                     (x) float32 -2763.22 -2760.68 -2758.14 -2755.6 ...
  * y                     (y) float32 -263.791 -261.251 -258.711 -256.172 ...
  * time                  (time) datetime64[ns] 2016-12-18 ...
    reftime               (time) datetime64[ns] ...
  * height_above_ground1  (height_above_ground1) float32 2.0
Attributes:
    long_name: Temperature @ Specified height level above ground
    units: K
    abbreviation: TMP
    grid_mapping: LambertConformal_Projection
    Grib_Variable_Id: VAR_0-0-0_L103
    Grib2_Parameter: [0 0 0]
    Grib2_Parameter_Discipline: Meteorological products
    Grib2_Parameter_Category: Temperature
    Grib2_Parameter_Name: Temperature
    Grib2_Level_Type: Specified height level above ground
    Grib2_Generating_Process_Type: Forecast

In [5]:

    

grid = nc[ncvar.grid_mapping]
grid


    

    
Out[5]:





<xarray.DataArray 'LambertConformal_Projection' ()>
array(0)
Attributes:
    grid_mapping_name: lambert_conformal_conic
    latitude_of_projection_origin: 25.0
    longitude_of_central_meridian: 265.0
    standard_parallel: 25.0
    earth_radius: 6371229.0
    _CoordinateTransformType: Projection
    _CoordinateAxisTypes: GeoX GeoY

In [6]:

    

lon0 = grid.longitude_of_central_meridian
lat0 = grid.latitude_of_projection_origin
lat1 = grid.standard_parallel
earth_radius = grid.earth_radius


    

In [7]:

    

import cartopy
import cartopy.crs as ccrs


    

In [8]:

    

isub = 10


    

In [9]:

    

ncvar.x


    

    
Out[9]:





<xarray.DataArray 'x' (x: 2145)>
array([-2763.21704102, -2760.67724609, -2758.13769531, ...,  2676.82666016,
        2679.36621094,  2681.90600586], dtype=float32)
Coordinates:
  * x        (x) float32 -2763.22 -2760.68 -2758.14 -2755.6 -2753.06 ...
Attributes:
    standard_name: projection_x_coordinate
    units: km
    _CoordinateAxisType: GeoX

In [10]:

    

#cartopy wants meters, not km
x = ncvar.x[::isub].data*1000.
y = ncvar.y[::isub].data*1000.


    

In [11]:

    

#globe = ccrs.Globe(ellipse='WGS84') #default
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=grid.earth_radius)

crs = ccrs.LambertConformal(central_longitude=lon0, central_latitude=lat0, 
                            standard_parallels=(lat0,lat1), globe=globe)


    

In [12]:

    

# find the correct time dimension name
for d in ncvar.dims:
    if "time" in d: 
        timevar = d


    

In [13]:

    

istep = -1
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
mesh = ax.pcolormesh(x,y,ncvar[istep,0,::isub,::isub].data.squeeze()-273.15, transform=crs,zorder=0, vmin=0, vmax=40)
fig.colorbar(mesh)
ax.coastlines(resolution='10m',color='black',zorder=1)
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
plt.title(nc[timevar].data[istep]);


    

    




C:\Users\rsignell\AppData\Local\Continuum\Miniconda3\envs\ioos3\lib\site-packages\matplotlib\colors.py:581: RuntimeWarning: invalid value encountered in less
  cbook._putmask(xa, xa < 0.0, -1)






    

In [ ]: