I'm trying to create a mask following this question's answer, so that I can color areas inside my shape file only. I will post the code example here (which is pretty much the same with the code in my previous question, but with the solution code included).
A quick summary of my code: It's a data of (lon, lat, tem). I transform the lon and lat into ccrs.PlateCaree() coordinate system (which is in meter and not degree anymore), and put the transformed coordinates in xp, yp. The 3 variables xp, yp, and tem are then interpolated and then I use pmeshcolor() to draw a color map (the same color map in my previous question that I left the link above). Now I want to limit the color map to be inside the shape file only (which is the grey lines). Code is as below:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import pandas as pd
import numpy as np
from metpy.interpolate import interpolate_to_grid, remove_nan_observations
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
canada_east = -95
canada_west = -101.8
canada_north = 52.8
canada_south = 48.85
central_lon = (canada_east + canada_west)/2
central_lat = (canada_north + canada_south)/2
crs = ccrs.LambertConformal(central_longitude = central_lon, central_latitude = central_lat)
lat = np.array([49.8134 50.904 50.698 49.095 49.436 49.9607 49.9601 49.356 50.116
49.402 52.3472 50.411 49.24 49.876 49.591 49.905 49.498 49.088
49.118 50.5947 49.3776 49.148 49.1631 51.358 49.826 50.4324 49.96
49.68 49.875 50.829 51.572])
lon = np.array([-100.3721 -97.273 -99.068 -97.528 -100.308 -98.9054 -98.6367
-99.248 -96.434 -100.93 -101.1099 -100.893 -100.055 -99.909
-97.518 -99.354 -98.03 -99.325 -99.054 -98.0035 -100.5387
-100.491 -97.1454 -100.361 -96.776 -99.4392 -97.7463 -97.984
-95.92 -98.111 -100.488])
tem = np.array([-8.45 -4.026 -5.993 -3.68 -7.35 -7.421 -6.477 -8.03 -3.834
-13.04 -4.057 -8.79 -6.619 -10.89 -4.465 -8.41 -4.861 -9.93
-7.125 -4.424 -11.95 -9.56 -3.86 -7.17 -4.193 -7.653 -4.883
-5.631 -3.004 -4.738 -8.81])
xp, yp, _ = crs.transform_points(ccrs.PlateCarree(), lon, lat ).T
xp, yp, tem = remove_nan_observations(xp, yp, tem)
alt_x, alt_y, data = interpolate_to_grid( xp, yp, tem, minimum_neighbors=2, search_radius=240000, interp_type = 'barnes', hres = 1000)
# Create the figure and grid for subplots
fig = plt.figure(figsize=(17, 12))
# Main ax
ax = plt.subplot(111, projection=crs)
ax.set_extent([canada_west, canada_east, canada_south, canada_north], ccrs.PlateCarree())
# Read shape file again using geopandas and create a mask
cat_gdf = geopandas.read_file('MB_AGregion_Perim_South.shp')
cat_gdf = cat_gdf.to_crs(epsg = 4326)
mask = shapely.vectorized.contains(cat_gdf.dissolve().geometry.item(), alt_x, alt_y)
print(mask)
which gives
[[False False False ... False False False]
[False False False ... False False False]
[False False False ... False False False]
...
[False False False ... False False False]
[False False False ... False False False]
[False False False ... False False False]]
Process finished with exit code 0
My alt_x and alt_y is in meter and my shape file geometries are in degree. I believe this is the reason why the mask contains all false value (which means masking everything). However, I'm unable to find a solution to this problem. I would appreciate a lot if someone can help me with this problem.
