Search code examples
projectioncartopy

Can't re-project precipitation data from Stereographic to PlateCarree() using Cartopy


I am trying to plot precipitation data from the National Weather Service. However, the data is by default set to a stereographic projection. I'd like to plot in a PlateCarree projection but I am having some difficulties. When I try and use the PlateCarree projection in Cartopy, it plots the maps but will not overlay the precipitation data. I'm assuming this means that I am not properly re-projecting the data from stereographic to PlateCarree. Is there anything specific I need to do in order to re-project the data correctly?

Here is the code that works with the stereographic projection:

'''

=====================
NWS Precipitation Map
=====================

Plot a 1-day precipitation map using a netCDF file from the National Weather Service.

This opens the data directly in memory using the support in the netCDF library to open
from an existing memory buffer. In addition to CartoPy and Matplotlib, this uses
a custom colortable as well as MetPy's unit support.



"""
###############################
# Imports
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.plots import USCOUNTIES
from metpy.units import masked_array, units
from netCDF4 import Dataset
import pandas as pd


###############################
# Download the data from the National Weather Service.
dt = datetime.utcnow() - timedelta(days=1)  # This should always be available
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)

###############################
# Pull the needed information out of the netCDF file
prcpvar = nc.variables['observation']
data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('in')
#data = data * 0.0393
x = nc.variables['x'][:]
y = nc.variables['y'][:]
proj_var = nc.variables[prcpvar.grid_mapping]

#%%
###############################
# Set up the projection information within CartoPy
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)


###############################
# Create the figure and plot the data
# create figure and axes instances
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(111, projection=proj)
#ax.set_extent([-75,-85,35,39])

#draw coastlines, state and country boundaries, edge of map.

ax.coastlines(resolution='10m')
ax.add_feature(cfeature.BORDERS.with_scale('10m'), linewidth=1.5)
ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=2.0)
ax.add_feature(USCOUNTIES.with_scale('500k'), edgecolor='black')

# draw filled contours.
clevs = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
          6.0, 8.0, 10., 20.0]
# In future MetPy
# norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)

cmap_data = [
    "#04e9e7",  # 0.01 - 0.10 inches
    "#019ff4",  # 0.10 - 0.25 inches
    "#0300f4",  # 0.25 - 0.50 inches
    "#02fd02",  # 0.50 - 0.75 inches
    "#01c501",  # 0.75 - 1.00 inches
    "#008e00",  # 1.00 - 1.50 inches
    "#fdf802",  # 1.50 - 2.00 inches
    "#e5bc00",  # 2.00 - 2.50 inches
    "#fd9500",  # 2.50 - 3.00 inches
    "#fd0000",  # 3.00 - 4.00 inches
    "#d40000",  # 4.00 - 5.00 inches
    "#bc0000",  # 5.00 - 6.00 inches
    "#f800fd",  # 6.00 - 8.00 inches
    "#9854c6",  # 8.00 - 10.00 inches
    "#fdfdfd"   # 10.00+
]

cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)

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

# add colorbar.
cbar = plt.colorbar(cs, orientation='vertical')
cbar.set_label(data.units)

time =  nc.creation_time[4:6]+'/'+nc.creation_time[6:8]+'/'+nc.creation_time[0:4]+'  '\
+nc.creation_time[9:11] +':'+  nc.creation_time[11:13] + "  UTC"
print(time)

ax.set_title('24 hr Precipitation (in)' + '\n for period ending ' + time, fontsize = 16, fontweight = 'bold' )

'''

However, when I change the projection lines to PlateCarree I run into the issue I described above. Does anyone have any advice on how to re-prroject this data?

Thanks


Solution

  • Just add the transform parameter to the ax.contourf method.

    # Imports
    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.plots import USCOUNTIES
    from metpy.units import masked_array, units
    from netCDF4 import Dataset
    import pandas as pd
    
    
    ###############################
    # Download the data from the National Weather Service.
    dt = datetime.utcnow() - timedelta(days=1)  # This should always be available
    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)
    
    ###############################
    # Pull the needed information out of the netCDF file
    prcpvar = nc.variables['observation']
    data = masked_array(prcpvar[:], units(prcpvar.units.lower())).to('in')
    #data = data * 0.0393
    x = nc.variables['x'][:]
    y = nc.variables['y'][:]
    proj_var = nc.variables[prcpvar.grid_mapping]
    
    #%%
    ###############################
    # Set up the projection information within CartoPy
    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)
    
    ###############################
    # Create the figure and plot the data
    # create figure and axes instances
    fig = plt.figure(figsize=(15, 15))
    
    pc_proj = ccrs.PlateCarree()
    #ax = fig.add_subplot(111, projection=proj)
    ax = fig.add_subplot(111, projection=pc_proj)
    
    
    ax.set_extent([-128,-65,25,52])
    #draw coastlines, state and country boundaries, edge of map.
    
    ax.coastlines(resolution='10m')
    #ax.add_feature(cfeature.BORDERS.with_scale('10m'), linewidth=1.5)
    #ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=2.0)
    #ax.add_feature(USCOUNTIES.with_scale('500k'), edgecolor='black')
    
    # draw filled contours.
    clevs = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
              6.0, 8.0, 10., 20.0]
    # In future MetPy
    # norm, cmap = ctables.registry.get_with_boundaries('precipitation', clevs)
    
    cmap_data = [
        "#04e9e7",  # 0.01 - 0.10 inches
        "#019ff4",  # 0.10 - 0.25 inches
        "#0300f4",  # 0.25 - 0.50 inches
        "#02fd02",  # 0.50 - 0.75 inches
        "#01c501",  # 0.75 - 1.00 inches
        "#008e00",  # 1.00 - 1.50 inches
        "#fdf802",  # 1.50 - 2.00 inches
        "#e5bc00",  # 2.00 - 2.50 inches
        "#fd9500",  # 2.50 - 3.00 inches
        "#fd0000",  # 3.00 - 4.00 inches
        "#d40000",  # 4.00 - 5.00 inches
        "#bc0000",  # 5.00 - 6.00 inches
        "#f800fd",  # 6.00 - 8.00 inches
        "#9854c6",  # 8.00 - 10.00 inches
        "#fdfdfd"   # 10.00+
    ]
    
    cmap = mcolors.ListedColormap(cmap_data, 'precipitation')
    norm = mcolors.BoundaryNorm(clevs, cmap.N)
    
    #cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm)
    # add transform args
    cs = ax.contourf(x, y, data, clevs, alpha = 0.5, cmap=cmap, norm=norm,transform=proj)
    
    # add colorbar.
    cbar = plt.colorbar(cs, orientation='vertical')
    cbar.set_label(data.units)
    time =  nc.creation_time[4:6]+'/'+nc.creation_time[6:8]+'/'+nc.creation_time[0:4]+'  '\
    +nc.creation_time[9:11] +':'+  nc.creation_time[11:13] + "  UTC"
    print(time)
    
    ax.set_title('24 hr Precipitation (in)' + '\n for period ending ' + time, fontsize = 16, fontweight = 'bold' )
    plt.savefig("prec_usa_pc.png")
    #plt.show()
    
    

    Below is the output fig.enter image description here