/usr/share/doc/python-mpltoolkits.basemap-doc/examples/polarmaps.py is in python-mpltoolkits.basemap-doc 1.0.7+dfsg-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | from __future__ import print_function
# make plots of etopo bathymetry/topography data on
# various map projections, drawing coastlines, state and
# country boundaries, filling continents and drawing
# parallels/meridians
# illustrates special-case polar-centric projections.
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
import matplotlib.pyplot as plt
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
etopo = np.loadtxt('etopo20data.gz')
lons = np.loadtxt('etopo20lons.gz')
lats = np.loadtxt('etopo20lats.gz')
print('min/max etopo20 data:')
print(etopo.min(),etopo.max())
# these are the 4 polar projections
projs = ['laea','stere','aeqd','ortho'] # short names
# long names
projnames = ['Lambert Azimuthal Equal Area','Stereographic','Azimuthal Equidistant','Orthographic']
# loop over hemispheres, make a 4-panel plot for each hemisphere
# showing all four polar projections.
for hem in ['North','South']:
if hem == 'South':
lon_0 = -130.
lon_0_ortho = lon_0 - 180.
lat_0 = -90.
# Lambert Azimuth bounding lat must not extend into opposite hem.
bounding_lat = -0.01
elif hem == 'North':
lon_0 = 130.
lon_0_ortho = lon_0
lat_0 = 90.
# Lambert Azimuth bounding lat must not extend into opposite hem.
bounding_lat = 0.01
# loop over projections, one for each panel of the figure.
fig = plt.figure(figsize=(8,8))
npanel = 0
for proj,projname in zip(projs,projnames):
npanel = npanel + 1
if hem == 'South':
projection = 'sp'+proj
elif hem == 'North':
projection = 'np'+proj
# setup map projection
# centered on Australia (for SH) or US (for NH).
if proj == 'ortho':
m = Basemap(projection='ortho',
resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
else:
m = Basemap(boundinglat=bounding_lat,lon_0=lon_0,\
resolution='c',area_thresh=10000.,projection=projection,round=True)
# compute native map projection coordinates for lat/lon grid.
x,y = m(*np.meshgrid(lons,lats))
ax = fig.add_subplot(2,2,npanel)
# make filled contour plot.
cs = m.contourf(x,y,etopo,np.linspace(-7500,4500,41),cmap=cm.GMT_haxby)
# draw coastlines.
m.drawcoastlines()
# draw parallels and meridians.
m.drawparallels(np.arange(-80.,90,20.))
#labels = [l,r,t,b]
m.drawmeridians(np.arange(0.,340.,30.),labels=[1,1,1,1],fontsize=7)
# draw boundary around map region.
m.drawmapboundary()
# draw title.
plt.title(hem+' Polar '+projname,y=1.05,fontsize=12)
print('plotting '+hem+' Polar '+projname+' basemap ...')
plt.show()
|