网站建设图片/seo外包公司费用

观察到，做相反的事情，即在海洋上放置一个光栅，在大陆上覆盖一个面具，这很容易。只需使用map.fillcontinents()。所以这个解决方案的基本思想是修改fillcontinents函数，使其在海洋上形成多边形。在

步骤如下：创建一个大的圆形多边形，覆盖整个地球。在

为map.coastpolygons数组中的每个形状创建一个多边形。在

使用shapely及其difference方法将陆地多边形的形状从圆中剪开。在

在顶部添加剩余的多边形，它们具有海洋的形状，具有高zorder。在

代码：from mpl_toolkits.basemap import Basemap

import numpy as np

from scipy import interpolate

from shapely.geometry import Polygon

from descartes.patch import PolygonPatch

def my_circle_polygon( (x0, y0), r, resolution = 50 ):

circle = []

for theta in np.linspace(0,2*np.pi, resolution):

x = r * np.cos(theta) + x0

y = r * np.sin(theta) + y0

circle.append( (x,y) )

return Polygon( circle[:-1] )

def filloceans(the_map, color='0.8', ax=None):

# get current axes instance (if none specified).

if not ax:

ax = the_map._check_ax()

# creates a circle that covers the world

r = 0.5*(map.xmax - map.xmin) # + 50000 # adds a little bit of margin

x0 = 0.5*(map.xmax + map.xmin)

y0 = 0.5*(map.ymax + map.ymin)

oceans = my_circle_polygon( (x0, y0) , r, resolution = 100 )

# for each coastline polygon, gouge it out of the circle

for x,y in the_map.coastpolygons:

xa = np.array(x,np.float32)

ya = np.array(y,np.float32)

xy = np.array(zip(xa.tolist(),ya.tolist()))

continent = Polygon(xy)

## catches error when difference with lakes

try:

oceans = oceans.difference(continent)

except:

patch = PolygonPatch(continent, color="white", zorder =150)

ax.add_patch( patch )

for ocean in oceans:

sea_patch = PolygonPatch(ocean, color="blue", zorder =100)

ax.add_patch( sea_patch )

########### DATA

x = [3.395833, 31.233333, 28.045556, 45.35 ]

y = [6.453056, 30.05, -26.204444, 2.033333]

z = [0, 90, 180, 270]

# set up orthographic map projection

map = Basemap(projection='ortho', lat_0=0, lon_0=20, resolution='l')

## Plot the cities on the map

map.plot(x,y,".", latlon=1)

# create a interpolated mesh and set it on the map

interpol_func = interpolate.interp2d(x, y, z, kind='linear')

newx = np.linspace( min(x), max(x) )

newy = np.linspace( min(y), max(y) )

X,Y = np.meshgrid(newx, newy)

Z = interpol_func(newx, newy)

map.pcolormesh( X, Y, Z, latlon=1, zorder=3)

filloceans(map, color="blue")

您好：