自动将Matplotlib Basemap居中于数据

Question

我希望您能提供一种自动将基础地图绘制在我的坐标数据中心的解决方案。

我已经实现了自动居中，但是结果区域比实际数据使用的区域大得多。我希望绘图被限制在绘图坐标而不是从纬度/经度边界划定的区域内。

我正在使用John Cook的代码计算球面上两点之间的距离。

首次尝试：

这是我最初的脚本。它导致宽度和高度太小，无法容纳数据区域，并且中心纬度（lat0）偏南。

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import sys
import csv
import spheredistance as sd


print '\n'
if len(sys.argv) < 3:
    print >>sys.stderr,'Usage:',sys.argv[0],'<datafile> <#rows to skip>'
    sys.exit(1)
print '\n'

dataFile = sys.argv[1]
dataStream = open(dataFile, 'rb')
dataReader = csv.reader(dataStream,delimiter='\t')
numRows = sys.argv[2]

dataValues = []
dataLat = []
dataLon = []

print 'Plotting Data From: '+dataFile

dataReader.next()
for row in dataReader:
    dataValues.append(row[0])
    dataLat.append(float(row[1]))
    dataLon.append(float(row[2]))

# center and set extent of map
earthRadius = 6378100 #meters
factor = 1.00

lat0new = ((max(dataLat)-min(dataLat))/2)+min(dataLat)
lon0new = ((max(dataLon)-min(dataLon))/2)+min(dataLon)

mapH = sd.distance_on_unit_sphere(max(dataLat),lon0new,
            min(dataLat),lon0new)*earthRadius*factor

mapW = sd.distance_on_unit_sphere(lat0new,max(dataLon),
            lat0new,min(dataLon))*earthRadius*factor

# setup stereographic basemap.
# lat_ts is latitude of true scale.
# lon_0,lat_0 is central point.
m = Basemap(width=mapW,height=mapH,
            resolution='l',projection='stere',\
            lat_0=lat0new,lon_0=lon0new)

#m.shadedrelief()
m.drawcoastlines(linewidth=0.2)
m.fillcontinents(color='white', lake_color='aqua')

#plot data points (omitted due to ownership)
#x, y = m(dataLon,dataLat)
#m.scatter(x,y,2,marker='o',color='k')

# draw parallels and meridians.
m.drawparallels(np.arange(-80.,81.,20.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180.,181.,20.), labels=[0,0,0,1], fontsize=10)
m.drawmapboundary(fill_color='aqua')

plt.title("Example")
plt.show()

- ryanjdillon

---

你应该将你的答案发布为一个回答（你也可以接受自己的答案），并将你的第二个问题作为一个新问题提出。 - bmu

只是重新审视这个。听起来不错。谢谢 bmu。 - ryanjdillon

Answer 1

生成了一些随机数据后，很明显我选择的边界与此投影不符（红线）。使用map.drawgreatcircle()，首先在随机数据投影缩放时可视化出我想要的边界。 我通过在最南端的纬度处使用经度差值来更正经度（蓝色水平线）。

我使用勾股定理解决垂直距离问题，已知最北端经度边界与中央最南端点之间的距离以及距离求解纬度范围（蓝色三角形）。

def centerMap(lats,lons,scale):
    #Assumes -90 < Lat < 90 and -180 < Lon < 180, and
    # latitude and logitude are in decimal degrees
    earthRadius = 6378100.0 #earth's radius in meters

    northLat = max(lats)
    southLat = min(lats)
    westLon = max(lons)
    eastLon = min(lons)

    # average between max and min longitude 
    lon0 = ((westLon-eastLon)/2.0)+eastLon

    # a = the height of the map
    b = sd.spheredist(northLat,westLon,northLat,eastLon)*earthRadius/2
    c = sd.spheredist(northLat,westLon,southLat,lon0)*earthRadius

    # use pythagorean theorom to determine height of plot
    mapH = pow(pow(c,2)-pow(b,2),1./2)
    arcCenter = (mapH/2)/earthRadius

    lat0 = sd.secondlat(southLat,arcCenter)

    # distance between max E and W longitude at most souther latitude
    mapW = sd.spheredist(southLat,westLon,southLat,eastLon)*earthRadius

    return lat0,lon0,mapW*scale,mapH*scale

lat0center,lon0center,mapWidth,mapHeight = centerMap(dataLat,dataLon,1.1)

在这种情况下，lat0（或纬度中心）是三角形高度的一半处的点，我使用John Cook的方法解决了这个问题，但是要解决一个未知坐标，同时知道第一个坐标（南边界的中位经度）和弧长（总高度的一半）。

def secondlat(lat1, arc):
    degrees_to_radians = math.pi/180.0
    lat2 = (arc-((90-lat1)*degrees_to_radians))*(1./degrees_to_radians)+90
    return lat2

更新: 可以使用pyproj Geod类方法geod.fwd()和geod.inv()，来实现上述函数以及两个坐标之间的距离计算，更高的精度。我在Erik Westra的《Python for Geospatial Development》中找到了这个信息，这是一个很好的资源。

更新: 我已经验证这也适用于Lambert Conformal Conic (lcc)投影。

- ryanjdillon

---

很棒的答案！基本上，你所做的是计算适合你数据坐标的最小直立球形矩形。请注意，你的代码目前只处理了北半球，尽管将其推广很简单。此外，可以查看geographiclib获取纯Python实现的大地测量计算。 - letmaik