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
1 回答
感谢您的回复，我确实通过使用我所在区域的edge-gridcells的坐标制作多边形来解决它 .
{
}