Orthographic projection Python

2019-02-20 20:14发布

问题:

I use orthographic projection to plot maps. I use this programm:

from mpl_toolkits.basemap import Basemap
import numpy as np    
import matplotlib.pyplot as plt   
import os, sys    
from sys import argv   
import pylab    
from mpl_toolkits.basemap import Basemap, shiftgrid    
from matplotlib import mpl     
from matplotlib import rcParams    
import matplotlib.pyplot as plt    
import matplotlib.mlab as mlab    
import matplotlib.patches as patches    
import matplotlib.path as path    
import matplotlib.dates as dt    
from numpy import linalg    
import netCDF4    
import time    
import datetime as d   
import sys    
import math    
from mpl_toolkits.axes_grid1 import make_axes_locatable   
from pylab import *


nc = netCDF4.Dataset ('tt.nc')    
latvar = nc.variables['lat']    
lat = latvar[:]    
lon = nc.variables['lon'][:]    
lat_0=30;lon_0=-25    
m1 = Basemap(projection='ortho',lon_0=-25,lat_0=30,resolution='l')    
m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution='l',\
    llcrnrx=0.,llcrnry=0.,urcrnrx=m1.urcrnrx/2.,urcrnry=m1.urcrnry/2.)

X, Y = m(lon, lat)    
O_x_1=nc.variables['O3']   
h=9    
lev=0    
minOzone=0    
maxOzone=40    
plotOzone = m.pcolor(X,Y,O_x_1[h,lev,:,:],vmin=minOzone,vmax=maxOzone)
ax=colorbar(plotOzone, shrink=0.8,norm=(0,40))    
m.drawcoastlines()   
m.drawparallels(np.arange(-90.,120.,30.))    
m.drawmeridians(np.arange(0.,420.,60.))    
plt.show()

What do I have to do to center my map on Europe ?

I have already played with lat_0 and lon_0 but that doesn't give what I want...

I can't add figures to show what I obtained and what I would like...

Thank you!

回答1:

The lat_0 and lon_0 are for setting the origin of the projection, not the extent of the map. Usually the place with the least distortions so you dont want the origin to deviate too much from the center of your area of interest. Basemap will automatically center the map around the origin if you dont specify an extent.

Centering the map (different from its origin) can be done if you know what extent (or boundingbox) you want to use. If you know the corner coordinates in your 'ortho' projection you could use the keywords from your example (llcrnrx etc.). I have had no luck with the 'llcrnrlon' keywords in Basemap 1.0.6, they seem to suggest that you can input the coordinates of your extent in geographic (lat/lon).

An alternative is to grab the axes and manually set an x- and y limit. A benefit is that you can do it after declaring the Basemap object which you can then use for coordinate transformation. An example:

from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(5,5))

m = Basemap(projection='ortho',lon_0=5,lat_0=35,resolution='l')

m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,15.))
m.drawmeridians(np.arange(0.,420.,30.))

# your extent in lat/lon (dec degrees)
ulx = -10
uly = 65
lrx = 65
lry = 35

# transform coordinates to map projection
xmin, ymin = m(ulx, lry)
xmax, ymax = m(lrx, uly)

# set the axes limits
ax = plt.gca()
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

Make sure the projection in the map declaration suites your needs, i have just picked an origin which falls within Europe.