Автоматически центрировать базовую карту matplotlib по данным

Мне нужно решение для автоматического центрирования графика базовой карты на моих координатных данных.

У меня есть вещи для автоматического центрирования, но результирующая область намного больше, чем область, фактически используемая моими данными. Я бы хотел, чтобы график был ограничен координатами графика, а не областью, нарисованной по границам широты и долготы.

Я использую код Джона Кука для вычисления расстояния между двумя точками на (предполагаемой идеальной) сфере .

Первая попытка

Вот сценарий, с которого я начал. Это приводило к тому, что ширина и высота были слишком малы для области данных, а центральная широта (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()

введите здесь описание изображения


person ryanjdillon    schedule 05.11.2012    source источник
comment
вы должны опубликовать свой ответ как ответ (вы также можете принять свой собственный ответ) и задать свой второй вопрос как новый вопрос.   -  person bmu    schedule 13.11.2012
comment
просто пересматриваю это. это звучит как план. спасибо бму.   -  person ryanjdillon    schedule 21.11.2012


Ответы (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 (или центр широты) в данном случае является точкой на полпути к высоте этого треугольника, которую я решил, используя метод Джона Кука, но для решения неизвестной координаты, зная первую координату (среднюю долготу в точке). южная граница) и длина дуги (половина общей высоты).

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(). Я нашел это в превосходном ресурсе Python for Geospatial Development Эрика Вестры.

Обновление: я убедился, что это также работает для проекций Ламберта, конформных конических (lcc).

person ryanjdillon    schedule 21.11.2012
comment
Отличный ответ! По сути, вы вычисляете наименьший вертикальный сферический прямоугольник, который соответствует вашим координатам данных. Обратите внимание, что ваш код на данный момент обрабатывает только северное полушарие, хотя его легко обобщить. Кроме того, взгляните на geographiclib для реализации геодезических вычислений на чистом Python. - person letmaik; 28.06.2014