Plot boundaries of specific region in basemap

1k Views Asked by At

I have defined a region of interest where I have tracked e.g. evaporation in time. Now i want to depict this region on a basemap plot by plotting only its boundaries. The region is defined as an (almost) global lat/lon array filled only with 1 at the Region's coordinates (like a land/sea mask, but for my specific region).

If people want to plot boundaries of a certain geometry they often refer to shapefiles (which i am unfamiliar with), but it seems an easy way to create a polygon and plot this polygon on a basemap. However, i cannot find info on creating a shapefile from an array similar to my 'Region array'.

What is your suggestion?

1

There are 1 best solutions below

1
On

Thanks for the responses, i indeed solved it with making a polygon with the coordinates of the edge-gridcells of my region.

{

import numpy as np
from netCDF4 import Dataset

def getRegion(latnrs,lonnrs, latitude, longitude, lsm):   
lsm_globe = lsm
for lat in range(0,len(latitude)):
    for lon in range(0,len(longitude)):
        if longitude[lon] < 1.5:
            lsm_globe[lat,lon] = 0.
        if longitude[lon] > 15:
            lsm_globe[lat,lon] = 0.
        if latitude[lat] < 48:
            lsm_globe[lat,lon] = 0.
        if latitude[lat] > 54:
            lsm_globe[lat,lon] = 0.

Region = lsm_globe



import matplotlib.path as mpath

coord_region = np.argwhere(Region>0)

lats = np.zeros(len(coord_region))
lons = np.zeros(len(coord_region))
for i in range(len(coord_region)):

    lats[i] = coord_region[i][0]
    lons[i] = coord_region[i][1]

uppergp = []
lowergp = []
for i in range(len(coord_region)-1):
    if lats[i] < lats[i+1]:
        uppergp.append( [lats[i], lons[i]] )
        lowergp.append( [lats[i+1], lons[i+1]] )
uppergp.append( [lats[-1], lons[-1]] )
lowergp.insert(0, [lats[0], lons[0]] )
lowergp.reverse()
boundgp = uppergp + lowergp


vertlist = []       
for i in range(len(boundgp)):
    vertlist.append( (longitude[int(boundgp[i][1])]+1.125/2., latitude[int(boundgp[i][0])]-1.125/2.))

verts = vertlist
# adding last vert to list to close poly
verts.append(verts[-1])


Path = mpath.Path
lineto = Path.LINETO
codes = [Path.MOVETO, Path.CLOSEPOLY]
for i in range(len(boundgp)-1):
    codes.insert(1, lineto)


boundgpcoord = mpath.Path(verts, codes)
return boundgpcoord, Region

}