[python #ULF-679183]: MetPy Mapping Question

Greetings!

This was a fun one to figure out. Here's code that worked for me:

    %matplotlib inline

    import cartopy.feature as cfeature
    import cartopy.crs as ccrs
    import matplotlib.path as mpath
    import matplotlib.pyplot as plt
    import metpy
    import numpy as np
    import xarray as xr

    ds = 
xr.open_dataset('/Users/rmay/Downloads/total_rainfall_20200101_20201231.nc',decode_times=False)
    data_var = ds.metpy.parse_cf('Total_precipitation')

    x = data_var.x
    y = data_var.y
    im_data = data_var.isel(time=0)

    # Get the states feature
    states = cfeature.STATES.with_scale('50m')

    # Find OK by seeing what intersects with a small box--can be replaced
    # if there's an easy OK shapefile/set of vertices available
    ok = next(feat.intersecting_geometries((-98, -97, 35, 36)))

    # Take the coordinates of that feature and turn into a Matplotlib Path
    clip_region = mpath.Path(ok.exterior.coords)

    fig = plt.figure(figsize=(14, 14), dpi=300)
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

    # Set the boundary of the map using that Path
    ax.set_boundary(clip_region, transform=ccrs.PlateCarree())

    # Plot image and specify it in its native projection
    # Control detail of interpolation with regrid_shape--which also controls 
speed
    ax.imshow(im_data, extent=(x.min(), x.max(), y.min(), y.max()),
              cmap='Greys', origin='lower' if y[0] < y[-1] else 'upper',
              transform=data_var.metpy.cartopy_crs, regrid_shape=2000)
    ax.add_feature(states, edgecolor='tab:red')
    ax.set_extent((-103.25, -94.25, 33.6, 37.1), crs=ccrs.PlateCarree())

You'll want to make sure you have pykdtree installed in your environment for 
better speed on the image remapping in imshow. Given the use for 3D printing, 
you *might* try `regrid_shape=4000`. I'll say that `regrid_shape=2000` took ~40 
seconds on my new-ish (non-M1) macbook pro.

Hope this helps!

Ryan

> Ryan,
> 
> 
> 
> I had a question about MetPy/Matplotlib plotting that I didn’t know if you
> had an answer or could point me possibly to someone that could. I have
> NetCDF files of radar estimated rainfall data in and around Oklahoma. I am
> trying to get a greyscale image of just the Oklahoma data so that I can
> import that into 3D printer software. Here is some simple code that I have
> been trying.
> 
> 
> 
> %matplotlib inline
> 
> import cartopy.feature as cfeature
> 
> import cartopy.crs as ccrs
> 
> import matplotlib.pyplot as plt
> 
> import xarray as xr
> 
> import metpy
> 
> 
> 
> ds = xr.open_dataset('total_rainfall_20200101_20201231.nc',
> decode_times=False)
> 
> data_var = ds.metpy.parse_cf('Total_precipitation')
> 
> 
> 
> x = data_var.x
> 
> y = data_var.y
> 
> im_data = data_var.isel(time=0)
> 
> 
> 
> fig = plt.figure(figsize=(14, 14), dpi=300)
> 
> ax = fig.add_subplot(1, 1, 1, projection=data_var.metpy.cartopy_crs)
> 
> #ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
> 
> 
> 
> ax.imshow(im_data, extent=(x.min(), x.max(), y.min(), y.max()),
> 
> cmap='Greys', origin='lower' if y[0] < y[-1] else 'upper')
> 
> ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='tab:red')
> 
> 
> 
> plt.show()
> 
> fig.savefig('_test_fig.png')
> 
> 
> 
> 
> 
> Do you know of anyway to mask or cut out just Oklahoma? Alternatively, can
> neighboring states just be individually painted white? Is there a way to
> reproject the map to Plate Carree to give the northern border of Oklahoma a
> straight line? Thank you so much for any advice or help.


Ticket Details
===================
Ticket ID: ULF-679183
Department: Support Python
Priority: Low
Status: Closed
===================
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata 
inquiry tracking system and then made publicly available through the web.  If 
you do not want to have your interactions made available in this way, you must 
let us know in each email you send to us.