Warum stimmen die kartografischen Projektionskoordinaten der südpolaren Basemap-Karte nicht mit denen von Datensätzen in derselben Projektion überein?

Question

Einige satellitengestützte Erdbeobachtungsprodukte liefern Breiten-/Längeninformationen, während andere die X/Y-Koordinaten innerhalb einer gegebenen Gitterprojektion liefern (und es gibt auch einige, die beides haben, siehe Beispiel). Mein Ansatz im zweiten Fall besteht darin, eine Grundkarten-Map zu erstellen, die dieselben Parameter (Projektion, Ellipsoid, Kartenursprung) wie vom Datenprovider hat, so dass die angegebenen X/Y-Werte den Grundkartenkoordinaten entsprechen. Wenn dies jedoch der Fall ist, stimmt die Geolokalisierung nicht mit anderen Datensätzen überein, einschließlich der Grundkartenküstenlinie. Ich habe dies mit drei verschiedenen Datensätzen aus verschiedenen vertrauenswürdigen Quellen erfahren. Für das minimale Beispiel verwende ich Landsat-Daten, die von der U.S. Geological Survey bereitgestellt werden, die sowohl X/Y-Koordinaten eines Südpolar-Stereografischen Gitters als auch die entsprechenden Lat/Lon-Koordinaten für alle vier Ecken des Bildes enthält.

Von einem Landsat Metafile wir bekommen (ID: LC82171052016079LGN00):

CORNER_UL_LAT_PRODUCT = -66,61490 CORNER_UL_LON_PRODUCT = -61,31816 CORNER_UR_LAT_PRODUCT = -68,74325 CORNER_UR_LON_PRODUCT = -58,04533 CORNER_LL_LAT_PRODUCT = -67,68721 CORNER_LL_LON_PRODUCT = -67,01109 CORNER_LR_LAT_PRODUCT = -69,94052 CORNER_LR_LON_PRODUCT = -64,18581 CORNER_UL_PROJECTION_X_PRODUCT = -2259300,000 CORNER_UL_PROJECTION_Y_PRODUCT = 1.236.000,000 CORNER_UR_PROJECTION_X_PRODUCT = -1.981.500,000 CORNER _UR_PROJECTION_Y_PRODUCT = 1236000.000 CORNER_LL_PROJECTION_X_PRODUCT = -2259300,000 CORNER_LL_PROJECTION_Y_PRODUCT = 958500,000 CORNER_LR_PROJECTION_X_PRODUCT = -1.981.500,000 CORNER_LR_PROJECTION_Y_PRODUCT = 958500,000

...

GROUP = PROJECTION_PARAMETERS MAP_PROJECTION = "PS" BEZUGSPUNKT = "WGS84" ELLIPSOID = "WGS84" VERTICAL_LON_FROM_POLE = 0,00000 TRUE_SCALE_LAT = -71,00000 FALSE_EASTING = 0 FALSCH_NORTHING = 0 GRID_CELL_SIZE_PANCHROMATIC = 15,00 GRI D_CELL_SIZE_REFLECTIVE = 30.00 GRID_CELL_SIZE_THERMAL = 30.00 orientation = "NORTH_UP" RESAMPLING_OPTION = "CUBIC_CONVOLUTION" END_GROUP = PROJECTION_PARAMETERS

von Basemap mit der rechten Kartenprojektion verwenden, sollten wir in der Lage sein, um die Ecke lat/lon Werte aus dem abzuleiten X/Y-Werte:

import numpy as np 
from mpl_toolkits.basemap import Basemap 

m=Basemap(resolution='h',projection='spstere', ellps='WGS84', boundinglat=-60,lon_0=180, lat_ts=-71) 

x_crn=np.array([-2259300,-1981500,-2259300,-1981500])# upper left, upper right, lower left, lower right 
y_crn=np.array([1236000, 1236000, 958500, 958500])# upper left, upper right, lower left, lower right 

x0, y0= m(0, -90) 
#Basemap coordinates at the south pole 
#note that (0,0) of the Basemap is in a corner of the map, 
#while other data sets use the south pole. 
#This is easy to take into account: 
lon_crn, lat_crn = m(x0-x_crn, y0-y_crn, inverse=True) 


print 'lon_crn: '+str(lon_crn) 
print 'lat_crn: '+str(lat_crn) 

Welche zurück:

lon_crn: [-61.31816102 -58.04532791 -67.01108782 -64.1858106 ] 
lat_crn: [-67.23548626 -69.3099076 -68.28071626 -70.47651326] 

Wie Sie sehen können, stimmen die Längengrade mit der gegebenen Genauigkeit mit denen aus der Metadatei überein, aber die Breitengrade sind zu niedrig.

Ich kann die Breiten nähern von:

lat_crn=(lat_crn+90.)*1.0275-90. 

Aber das ist nicht wirklich befriedigend.

Dies ist, wie das Bild befindet, wenn unter Verwendung des X/Y Ecke von der Metadatei-Koordinaten (in rot die Basiskarte drawcoastlines()): und dies ist, wie es aussieht wie die Ecke unter Verwendung von lat/lon:

In diesem Fall kann ich einfach die Lat/Lon-Koordinaten verwenden, aber wie bereits erwähnt, gibt es Datensätze (wie this), die nur durch X/Y-Koordinaten zur Verfügung gestellt werden. Daher ist es sehr wichtig, sich auf die Basemap-Projektion zu verlassen. Ich weiß, dass es andere Module gibt, um die Daten als mögliche Umgehungslösung neu zu projizieren, aber es sollte ohne andere Module funktionieren und eine Neuprojektion könnte selbst Fehler verursachen.

Da dieses Problem bei verschiedenen Datensätzen auftritt, glaube ich gerne, dass es sich um einen Fehler im Basemap-Modul handelt, aber ich könnte auch immer wieder denselben Fehler machen oder falsche Erwartungen haben.

Answer 1

Ich habe etwas experimentiert und es scheint, als ob die Änderung lat_ts hat keinen Effekt mit projection='spstere'. In der Tat scheint es, als ob der Projektionsspielraum implizit als lat_ts=-90. angenommen wird, unabhängig davon, welchen Wert Sie zuweisen.

hatte ich mehr Erfolg projection='stere' stattdessen verwenden, so dass Sie die Basiskarte in Ihrem Beispiel konstruieren würde, wie folgt:

m=Basemap(width=5400000., height=5400000., projection='stere', 
      ellps='WGS84', lon_0=180., lat_0=-90., lat_ts=-71.) 

Sie bevorzugen anstelle der Breite und Höhe der Breite und Länge der Ecken einstellen der Handlung für Ihre Anwendung.