Search code examples
pythonmatplotlibcoordinatesprojection

Matplotlib: Focus on specific lon/lat using spstere projection

I am trying to focus my map on a specific area of Antarctica using 'spstere' projection from the matplotlib package in Python. I am able to plot the whole of Antarctica but this time I want to 'zoom' in and have a closer look at a specific area of the continent.

Similar examples using other projections (Pyplot contour plot - clabel spacing; http://matplotlib.org/basemap/api/basemap_api.html; https://matplotlib.org/basemap/users/examples.html) are available online but I have not been able to apply those to the 'spstere' projection over Antarctica.

I basically want to focus my map on the region of the Antarctic Peninsula, which spans roughly from

llcrnrlon=-100,urcrnrlon=-30,llcrnrlat=-90,urcrnrlat=-55.0

I have tried to use this code with the 'spstere' proj but python only takes into account boundinglat and lon_0. I've tried to change the values for boundinglat and lon_0 but it does not work either.

Any idea how I could go about? I have also tried using other projections such as 'cyl' but instead of getting a nice square like the 'spstere' proj, I get a horizontal rectangle.

m = Basemap(projection='cyl',lon_0=0,lat_0=0,\
      llcrnrlon=-180,urcrnrlon=180,llcrnrlat=-90,urcrnrlat=-55.0,resolution='c')

Any help would be highly appreciated!

Solution

  • Using the Polar Stereographic Projection 'spstere', you can get the antarctic region by using e.g. boundinglat=-60:

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

m = Basemap(projection='spstere',boundinglat=-60,lon_0=180,resolution='c')
m.drawcoastlines()

plt.show()

    enter image description here

    Note that 'spstere' is always centered at the south pole.
    In order to have a map, which is not centered at the south pole, you need to use the "stere" projection. Setting the corners for the "stere" projection is not straigt forward.

    One may therefore use a plot in 'spstere' projection and find some points which would enclose the region of interest. In this case e.g.

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

m = Basemap(projection='spstere',boundinglat=-50,
            lon_0=180+(-100+-30)/2.,resolution='c')

m.drawmeridians(np.arange(0,360,30),labels=[1,1,1,0])
m.drawparallels(np.arange(-90,90,5))
m.drawcoastlines()

xll, yll = m(-150,-70) # <-- find those points by looking at meridians and parallels
xur, yur = m(-30,-55)
m.scatter([xll,xur], [yll, yur], c="crimson")
plt.show()

    enter image description here

    Using those points, (-150,-70, -30,-55), as the corners of the map, you can then plot a map using 'stere' projection.

    m = Basemap(projection='stere',resolution='c',
            lat_0=-90, lon_0=(-100+-30)/2.,  lat_ts=(-90.+-55.)/2.,
            llcrnrlon=-150,urcrnrlon=-30,llcrnrlat=-70,urcrnrlat=-55)

    enter image description here

    If this heuristic method is not wanted, you may automate this procedure by creating a dummy map in 'spstere' projection, calculate the coordinates from the rectangle in question (llcrnrlon=-100,urcrnrlon=-30,llcrnrlat=-90,urcrnrlat=-55.0) and create a new basemap in stere projection with them. The function below is taken from the ActiveState site (author PG).

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

def polar_stere(lon_w, lon_e, lat_s, lat_n, **kwargs):
    '''Returns a Basemap object (NPS/SPS) focused in a region.
    lon_w, lon_e, lat_s, lat_n -- Graphic limits in geographical coordinates.
                                  W and S directions are negative.
    **kwargs -- Aditional arguments for Basemap object.
    '''
    lon_0 = lon_w + (lon_e - lon_w) / 2.
    ref = lat_s if abs(lat_s) > abs(lat_n) else lat_n
    lat_0 = np.copysign(90., ref)
    proj = 'npstere' if lat_0 > 0 else 'spstere'
    prj = Basemap(projection=proj, lon_0=lon_0, lat_0=lat_0,
                          boundinglat=0, resolution='c')
    lons = [lon_w, lon_e, lon_w, lon_e, lon_0, lon_0]
    lats = [lat_s, lat_s, lat_n, lat_n, lat_s, lat_n]
    x, y = prj(lons, lats)
    ll_lon, ll_lat = prj(min(x), min(y), inverse=True)
    ur_lon, ur_lat = prj(max(x), max(y), inverse=True)
    return Basemap(projection='stere', lat_0=lat_0, lon_0=lon_0,
                           llcrnrlon=ll_lon, llcrnrlat=ll_lat,
                           urcrnrlon=ur_lon, urcrnrlat=ur_lat, **kwargs)

llcrnrlon=-100
urcrnrlon=-30
llcrnrlat=-90
urcrnrlat=-55.0
m = polar_stere(llcrnrlon, urcrnrlon, llcrnrlat, urcrnrlat)

m.drawmeridians(np.arange(0,360,30),labels=[1,1,1,0])
m.drawparallels(np.arange(-90,90,30),labels=[1,1,1,1])
m.drawcoastlines()

plt.show()

    enter image description here