Search code examples
pythonnumpycartopy

How do I label a section of a cartopy map as NaN?

I have attached the map I am working with Map1. I am looking to change the values in the large section to the right of the blue line (labelled in the second image I have provided LabelledMap1) from int to NaN. Is there a way to choose a specific lat/lon coordinate and NaN all the values surrounding it? Or is there another method that would be more effective? A colleague recommended a flexible barrier to block out that section of the map, but I am not sure how to go about coding that.

How my map is produced:

fig = plt.figure(figsize=(10,10)) ax = fig.add_subplot(111,projection=rot) plot_cartopy(ax,lon,lat) plt.title("Month = May")

# set up the cartopy map

plt.contourf(lon,
            lat,
            new[4,:,:], 
            transform=ccrs.PlateCarree(), # tell matplotlib this is a map 
            cmap=cmocean.cm.thermal, # colormap 
            zorder=4 # zorder needs to be set so that it knows to plot it on top of the cartopy map  ) plt.clim(0,100) plt.colorbar()

Where "new" is a 3D array with the dimensions [month, x-coordinates, y-coordinates].

Thank you!

Solution

  • The following code gives you some examples of how to change values to NaN using numpy.where function.

    Since I don't have your data and don't know your map extents or line function, I randomly generate my data and line. I also generate an extent that can reproduce better your map extent, but the code should be valid for any map extent you have. For you to better understand there are plotting examples for both global map and subset extent.

    In the last example I change to NaN the values to the Right of a custom line.

    import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import cartopy
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point

# Create array
lat = np.linspace(-90,90,80)
lon = np.linspace(-180,180,120)
new = np.random.normal(size=[12,len(lat),len(lon)])

# Set map extent because my example data is global
# you wouldn't need to do this cause your data has already the proper extent
extent=[-20, 15, 5, 25] # (20°W to 15°E, 5°N to 25°N)

# Example global plot and extent (for new[4,...] as in your case)
def plot(data,title=""):
    fig = plt.figure()
    axes = [fig.add_subplot(2,1,i+1,projection=ccrs.PlateCarree()) for i in range(2)]
    # Get cyclic lon data
    plottedData,cyclicLon = add_cyclic_point(data, coord=lon, axis=-1)
    # Set extent and create rectangle box
    rect = patches.Rectangle((extent[0], extent[2]), extent[1]-extent[0], extent[3]-extent[2], 
        linewidth=1, edgecolor='r', facecolor='none', linestyle="--",zorder=3)
    for i,ax in enumerate(axes):
        ax.add_feature(cartopy.feature.COASTLINE)
        ax.add_feature(cartopy.feature.LAND, facecolor=[0.7,0.7,0.7,1])
        ax.add_feature(cartopy.feature.OCEAN, facecolor='white',zorder=2)
        im=ax.contourf(cyclicLon,lat,plottedData[4,...], cmap=plt.cm.Spectral, levels=20)
        gl=ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, alpha=0.5, linestyle='--')
        gl.right_labels = False
        gl.top_labels = False
        plt.colorbar(im)
    # Set extent to first ax and add rectangle box to second ax
    axes[0].add_patch(rect)
    axes[1].set_extent(extent)
    fig.suptitle(title)
    return axes

# Set spatial conditions on an array with ndim(t,lat,lon)
def latCondition(cond):
    return np.broadcast_to(cond,new.transpose(0,2,1).shape).transpose(0,2,1)
def lonCondition(cond):
    return np.broadcast_to(cond,new.shape)
def spatialCondition(latCond=True,lonCond=True):
    return np.logical_and(latCondition(latCond),lonCondition(lonCond))

# Base map
plot(new,title="Base Map")

####
# CHANGE SOME VALUES TO NaN USING numpy.where
####

# EX.1 Change values based on specific value
# Change points whose value > specificValue
specificValue=1
changed=np.where(new>specificValue,np.nan,new)
plot(changed,title=f"Values > {specificValue} changed to NaN")

# EX.2 Change values based on longitude:
# Change values whose longitude > longValue
longValue = -8
condition = spatialCondition(lonCond=lon > longValue)
changed=np.where(condition,np.nan,new)
plot(changed,title=f"Values with lon > {longValue} changed to NaN")

# EX.2 Change values based on longitude and latitude:
# Change values whose longitude is between longValue1 and longValue2 AND latitude < latValue
longValue1 = -115
longValue2 = 10
latValue = 12
condition = spatialCondition(latCond=lat<latValue,
    lonCond=np.logical_and(lon>longValue1, lon<longValue2))
changed=np.where(condition,np.nan,new)
plot(changed,title=f"Values with {longValue1} < lon < {longValue2} and\n lat < {latValue} changed to NaN")


# Change values based on custom line

# Generate custom line

# In my case to avoid having to interpolate the line over the map's grid,
# I generate the line over the same grid. 
# Also, to have something similar to what you have, 
# I generate my line as a vertical line, with each longitude value
# randomly chosen within the longitude extent that we selected.

lonExtent = lon[np.logical_and(lon >= extent[0],lon <= extent[1])]
lineLon = np.random.choice(lonExtent,size=len(lat))
line = [lineLon,lat]

# Get representation of custom line over the map
haxes=plot(new,title="Custom Line")
for hax in haxes:
    hax.plot(*line,
        color='blue', linewidth=2,
        transform=ccrs.PlateCarree())

# Change values to the RIGHT of the line to NaN
# So the condition is --> data's longitude >= lineLon[i] for i in range(lat)
condition = np.vstack([lon >= l for l in lineLon]) #2D Condition
changed=np.where(condition,np.nan,new)
haxes=plot(changed,title="Values to the RIGHT of the line changed to NaN")
for hax in haxes:
    hax.plot(*line,
        color='blue', linewidth=2,
        transform=ccrs.PlateCarree())
# Scatterplot to show NaN values
haxes[1].scatter(lon[np.where(np.isnan(changed[4,...]))[1]],lat[np.where(np.isnan(changed[4,...]))[0]],color='black',s=3,zorder=4)

    The reason why the results might seem not to be precise (for example the NaN values don't allign with the custom line in the plotting) is because of matplotlib contourf interpolation. Since matplotlib does not interpolate points if any of the values are NaNs, any point close to a NaN will result in that space not being interpolated for the visual interpolation. In the last plot I added a scatterplot of the NaN points to show you they are exactly the ones to the right of the custom line.

    Let me know if that helps.