Почему Matplotlib + Basemap показывает острова?

Question

Я работал с matplotlib и basemap, чтобы показать некоторую информацию о Нью-Йорке. До сих пор я следовал этому руководству , но столкнулся с проблемой. Я пытаюсь показать остров Манхэттен в своей визуализации, но не могу понять, почему базовая карта не отображает его как остров.

Вот визуализация, которую дает мне базовая карта:

Вот скриншот ограничивающей рамки, которую я использую:

А вот код, который генерирует изображение:

wl = -74.04006
sl = 40.683092
el = -73.834067
nl = 40.88378

m = Basemap(resolution='f', # c, l, i, h, f or None
            projection='merc',
            area_thresh=50,
            lat_0=(wl + sl)/2, lon_0=(el + nl)/2,
            llcrnrlon= wl, llcrnrlat= sl, urcrnrlon= el, urcrnrlat= nl)

m.drawmapboundary(fill_color='#46bcec')
m.fillcontinents(color='#f2f2f2',lake_color='#46bcec')
m.drawcoastlines()
m.drawrivers()

Я думал, что он мог бы рассмотреть воду между рекой, но m.drawrivers(), похоже, не мог это исправить. Очевидно, что любая помощь чрезвычайно ценится.

Заранее спасибо!

Answer 1

Один из способов получить более качественную базовую карту для ваших графиков - это построить ее из плиток веб-карты с соответствующим уровнем масштабирования. Здесь я продемонстрирую, как получить их с серверов openstreetmap веб-карт. В этом случае я использую уровень масштабирования 10 и объединяю 2 фрагмента карты в один массив изображений. Один из недостатков - экстент объединенного изображения всегда больше, чем запрашиваемые нами значения. Вот рабочий код:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import math
import urllib2
import StringIO
from PIL import Image

# === Begin block1 ===
# Credit: BerndGit, answered Feb 15 '15 at 19:47. And ...
# Source: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

def deg2num(lat_deg, lon_deg, zoom):
    '''Lon./lat. to tile numbers'''
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

def num2deg(xtile, ytile, zoom):
    '''Tile numbers to lon./lat.'''
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lat_deg, lon_deg)   # NW-corner of the tile. 

def getImageCluster(lat_deg, lon_deg, delta_lat, delta_long, zoom):
    # access map tiles from internet
    # no access/key or password is needed
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    # useful snippet: smurl.format(zoom, xtile, ytile) -> complete URL
    # x increases L-R; y Top-Bottom
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)  # get tile numbers (x,y)
    xmax, ymin =deg2num(lat_deg+delta_lat, lon_deg+delta_long, zoom)

    # PIL is used to build new image from tiles
    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl = smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = urllib2.urlopen(imgurl).read()
                # TODO: study, what these do?
                tile = Image.open(StringIO.StringIO(imgstr))
                Cluster.paste(tile, box=((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None
    return Cluster
# ===End Block1===

# Credit to myself
def getextents(latmin_deg, lonmin_deg, delta_lat, delta_long, zoom):
    '''Return LL and UR, each with (long,lat) of real extent of combined tiles.
       latmin_deg: bottom lat of extent
       lonmin_deg: left long of extent
       delta_lat: extent of lat
       delta_long: extent of long, all in degrees
    '''
    # Tile numbers(x,y): x increases L-R; y Top-Bottom
    xtile_LL, ytile_LL = deg2num(latmin_deg, lonmin_deg, zoom)  #get tile numbers as specified by (x, y)
    xtile_UR, ytile_UR = deg2num(latmin_deg + delta_lat, lonmin_deg + delta_long, zoom)
    # from tile numbers, we get NW corners
    lat_NW_LL, lon_NW_LL = num2deg(xtile_LL, ytile_LL, zoom)
    lat_NW_LLL, lon_NW_LLL = num2deg(xtile_LL, ytile_LL+1, zoom) # next down below

    lat_NW_UR, lon_NW_UR = num2deg(xtile_UR, ytile_UR, zoom)
    lat_NW_URR, lon_NW_URR = num2deg(xtile_UR+1, ytile_UR, zoom) # next to the right
    # get extents
    minLat = lat_NW_LLL
    minLon = lon_NW_LL
    maxLat = lat_NW_UR
    maxLon = lon_NW_URR
    return (minLon, maxLon, minLat, maxLat)  # (left, right, bottom, top) in degrees

# OP's values of extents for target area to plot
# some changes here (with larger zoom level) may lead to better final plot
wl = -74.04006
sl = 40.683092
el = -73.834067
nl = 40.88378

lat_deg = sl
lon_deg = wl
d_lat = nl - sl
d_long = el - wl
zoom = 10   # zoom level

# Acquire images. The combined images will be slightly larger that the extents
timg = getImageCluster(lat_deg, lon_deg, d_lat, d_long, zoom)

# This computes real extents of the combined tile images, and get (left, right, bottom, top)
latmin_deg, lonmin_deg, delta_lat, delta_long = sl, wl, nl-sl, el-wl
(left, right, bottom, top) = getextents(latmin_deg, lonmin_deg, delta_lat, delta_long, zoom)   #units: degrees

# Set Basemap with proper parameters
m = Basemap(resolution='h',  # h is nice
            projection='merc',
            area_thresh=50,
            lat_0=(bottom + top)/2, lon_0=(left + right)/2,
            llcrnrlon=left, llcrnrlat=bottom, urcrnrlon=right, urcrnrlat=top)

fig = plt.figure()
fig.set_size_inches(10, 12)
m.imshow(np.asarray(timg), extent=[left, right, bottom, top], origin='upper' )
m.drawcoastlines(color='gray', linewidth=3.0)  # intentionally thick line
#m.fillcontinents(color='#f2f2f2', lake_color='#46bcec', alpha=0.6)
plt.show()

Надеюсь, это поможет. Полученный участок:

Редактировать

Обрезать изображение, чтобы получить точную площадь для построения, не сложно. Модуль PIL может справиться с этим. Также работает срез массива Numpy.