我想要一个解决方案来自动将 basemap 绘图置于我的坐标数据的中心。
我有自动居中的东西，但结果区域比我的数据实际使用的区域大得多。 我希望地 block 以地 block 坐标为界，而不是从纬度/经度边界绘制的区域。
我正在使用 John Cook's code用于计算(假定的完美)球体上两点之间的距离。
第一次尝试
这是我开始使用的脚本。这导致数据区域的宽度和高度太小，中心纬度 (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()
请您参考如下方法:
生成一些随机数据后，很明显我选择的边界不适用于此投影(红线)。使用 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 Cooks 方法求解，但是在知道第一个坐标(中位数)的情况下求解未知坐标南界经度)和弧长(总高度的一半)。
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) 投影。
