Цветная карта matplotlib с использованием бикубической интерполяции

Я знаю, что matplotlib и scipy могут выполнять бикубическую интерполяцию: http://matplotlib.org/examples/pylab_examples/image_interp.html http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

Я также знаю, что можно нарисовать карту мира с помощью matplotlib: http://matplotlib.org/basemap/users/geography.html http://matplotlib.org/basemap/users/examples.html http://matplotlib.org/basemap/api/basemap_api.html

Но могу ли я сделать бикубическую интерполяцию на основе 4 точек данных и только окрасить массу земли?

Например, используя их для 4 точек данных (долгота и широта) и цветов:

Lagos: 6.453056, 3.395833; red HSV 0 100 100 (or z = 0)
Cairo: 30.05, 31.233333; green HSV 90 100 100 (or z = 90)
Johannesburg: -26.204444, 28.045556; cyan HSV 180 100 100 (or z = 180)
Mogadishu: 2.033333, 45.35; purple HSV 270 100 100 (or z = 270)

Я думаю, что нужно сделать бикубическую интерполяцию по широте широт и долгот, а затем добавить океаны, озера и реки поверх этого слоя? Я могу сделать это с помощью drawmapboundary. На самом деле для этого есть опция maskoceans: http://matplotlib.org/basemap/api/basemap_api.html#mpl_toolkits.basemap.maskoceans

Я могу интерполировать данные следующим образом:

xnew, ynew = np.mgrid[-1:1:70j, -1:1:70j]
tck = interpolate.bisplrep(x, y, z, s=0)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

Или с помощью scipy.interpolate.interp2d: http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

Здесь объясняется, как преобразовать в координаты проекции карты: http://matplotlib.org/basemap/users/mapcoords.html

Но мне нужно выяснить, как это сделать для расчетной поверхности вместо отдельных точек. На самом деле есть пример такой топографической карты с использованием внешних данных, которые я должен иметь возможность реплицировать: http://matplotlib.org/basemap/users/examples.html

P.S. Я не ищу полного решения. Я бы предпочел решить это сам. Скорее, я ищу предложения и подсказки. Я использую gnuplot более 10 лет и перешел на matplotlib в течение последних нескольких недель, поэтому, пожалуйста, не предполагайте, что я знаю даже самые простые вещи о matplotlib.

Ответ 1

Я думаю, что это то, что вы ищете (грубо). Обратите внимание, что важнейшие вещи маскируют массив данных до того, как вы построите pcolor и перейдете в параметр hsv colormap (Docs: cmap для pcolormesh и доступные цветовые карты).

Я сохранил код для построения карт, близких к примерам, поэтому его следует легко отслеживать. По той же причине я сохранил ваш код интерполяции. Обратите внимание, что интерполяция линейна, а не кубическая - kx=ky=1 - потому что вы не даете достаточного количества точек для кубической интерполяции (вам понадобится не менее 16 - scipy будет жаловаться на меньшее, говоря, что "m must be >= (kx+1)(ky+1)", хотя ограничение не упоминается в документации).

Я также расширил диапазон вашего meshgrid и сохранял в lat/lon для x и y повсюду.

Код

from mpl_toolkits.basemap import Basemap,maskoceans
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
# set up orthographic map projection with
# perspective of satellite looking down at 0N, 20W (Africa in main focus)
# use low resolution coastlines.
map = Basemap(projection='ortho',lat_0=0,lon_0=20,resolution='l')

# draw coastlines, country boundaries
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
# Optionally (commented line below) give the map a fill colour - e.g. a blue sea
#map.drawmapboundary(fill_color='aqua')
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))


data = {'Lagos': (6.453056, 3.395833,0),
        'Cairo': (30.05, 31.233333,90),
        'Johannesburg': (-26.204444, 28.045556,180),
        'Mogadishu': (2.033333, 45.35, 270)}

x,y,z = zip(*data.values())

xnew, ynew = np.mgrid[-30:60:0.1, -50:50:0.1]

tck = interpolate.bisplrep(x, y, z, s=0,kx=1,ky=1)
znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)
znew = maskoceans(xnew, ynew, znew)

col_plot = map.pcolormesh(xnew, ynew, znew, latlon=True, cmap='hsv')
plt.show()

Выход

Ответ 2

Соблюдайте, что делать наоборот, то есть наносить растр на море и накладывать маску на континенты, легко, как пирог. Просто используйте map.fillcontinents(). Таким образом, основная идея этого решения состоит в том, чтобы изменить функцию fillcontinents так, чтобы она закладывала многоугольники над океанами.

Шаги:

Создайте большой кругподобный многоугольник, покрывающий весь земной шар.
Создайте многоугольник для каждой фигуры в массиве map.coastpolygons.
Отрежьте форму многоугольника landmass от круга, используя 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")

Вуаля: