How to plot a shapefile centered in the Pacific wi

问题:

When plotting with Basemap's readshapefile, if the defined map is centered anywhere else than the longitudinal center of the shapefile, only a portion of it it's plotted. Here's an example using Natural Earth's coastlines:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

shpf = './NaturalEarth/ne_50m_land/ne_50m_land'

fig, ax = plt.subplots(nrows=1, ncols=1, dpi=100)

m = Basemap(
    ax = ax,
    projection = 'cyl',
    llcrnrlon = 0, llcrnrlat = -90,
    urcrnrlon = 360, urcrnrlat = 90
) 

m.readshapefile(shpf,'ne_50m_land')

m.drawmeridians(np.arange(0,360,45),labels=[True,False,False,True])

Which produces:

Is there a workaround for this with Basemap or Python? I know some people re-center the shapefile in QGIS or similar, but it seems unpractical to do so every time you create a new map, and my QGIS skills are extremely basic.

回答1:

One way to do it would be to tell readshapefile not to plot the coastlines directly and then to manipulate the line segments before plotting them yourself. Here an example based on your use case:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

shpf = 'shapefiles/ne_50m_land'

fig, ax = plt.subplots(nrows=1, ncols=1, dpi=100)

m = Basemap(
    ax = ax,
    projection = 'cyl',
    llcrnrlon = 0, llcrnrlat = -90,
    urcrnrlon = 360, urcrnrlat = 90
) 

m.readshapefile(shpf,'ne_50m_land', drawbounds = False)

boundary = 0.0

for info, shape in zip(m.ne_50m_land_info, m.ne_50m_land):
    lons, lats = map(np.array, zip(*shape))

    sep = (lons <= boundary).astype(int)
    roots = np.where(sep[:-1]+sep[1:] == 1)[0]+1
    lower = np.concatenate([[0],roots]).astype(int)
    upper = np.concatenate([roots,[len(lons)]]).astype(int)

    for low, high in zip(lower,upper):
        lo_patch = lons[low:high]
        la_patch = lats[low:high]
        lo_patch[lo_patch<0] += 360
        x,y = m(lo_patch,la_patch)
        ax.plot(x,y,'k',lw=0.5)

m.drawmeridians(np.arange(0,360,45),labels=[True,False,False,True])

plt.show()

In the example above, I iterate through the line segments of the shape file the way it is explained in the Basemap documentation. First I thought it would be enough to just add 360 to each point with a longitude smaller 0, but then you would get horizontal lines whenever a coast line crosses the 0 degree line. So, instead, one has to cut the lines into smaller segments whenever such a crossing appears. This is quite easily accomplished with numpy. I then use the plot command to draw the coast lines. If you want to do something more complex have a look at the Basemap documentation.

The final result looks like this:

Hope this helps.