ARTICLE AD BOX
I am trying to create a masking array for global map using shapefiles. I managed to get the right shapefile to work in a way that I can now plot an area that I want using matplotlib polygon. But, I just can't figure out how to map this shape into an array of my choice.
For example, I would like to map the shape of North Atlantic Ocean into the 0.1 degree resolution global grid (see code below). Is it possible to somehow get the coordinates of the Polygon that I have created and map it into the grid array in the code?
The shapefile is from https://marineregions.org/downloads.php "Global Oceans and Seas".
from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection from matplotlib.patches import PathPatch import numpy as np idir = "shp.files" ifn = "goas_v01" ifp = "{0:s}/{1:s}".format(idir,ifn) sres = 0.1 decim = int(-np.log10(sres)) lat_grid = np.arange( -90, 90+sres,sres).round(decim) lon_grid = np.arange(-180,180 ,sres).round(decim) lat_size = len(lat_grid) lon_size = len(lon_grid) grid = np.zeros((lat_size,lon_size)) mapargs = { "llcrnrlon" :-180 , "urcrnrlon" : 180 , "llcrnrlat" : -90 , "urcrnrlat" : 90 , "projection":"cyl", "lat_0" : 0 , "lon_0" : 0 , } fig,ax = plt.subplots(1,1) patches = [] map = Basemap(**mapargs) map.drawmapboundary() map.drawcoastlines() shpf = map.readshapefile("{0:s}".format(ifp),"{0:s}".format(ifn),drawbounds = False) infos = map.goas_v01_info shps = map.goas_v01 msize = len(infos) trg_nm = "North Atlantic Ocean" for m in range(msize): info = infos [m] rngnum = info ["RINGNUM"] nm = info ["name"] shp = shps [m] if (nm == trg_nm and rngnum == 1): patches.append(Polygon(np.array(shp), closed=True)) ax.add_collection(PatchCollection(patches, facecolor= 'm', edgecolor='k', linewidths=1., zorder=2)) plt.show()So, I would like to map the purple area (see picture)
into the grid array in dimensions of (1801,3600) as value 1 for each points.
