Geotransformation für polare Stereographie?

16

Ich arbeite derzeit daran, CANGRID-Klimadaten (bereitgestellt als Surfer Grid-ASCII-Dateien ".grd") in ArcGIS zu importieren. Das Raster ist 95 Zeilen mal 125 Spalten groß. Die Metadaten enthalten die Länge des Ursprungs (linke untere Ecke), die Zellengröße (50 km) sowie die Notenprojektion als polare Stereografie mit Mittelmeridian (110 Grad W) und den Breitengrad des Ursprungs (60 Grad N).

Nachdem ich zum ersten Mal erfolglos versucht habe, die GRD-Datei in Raster-Dateien wie .ascii und .flt zu konvertieren, konnte ich GDAL zum Festlegen der Ausdehnung und Projektion verwenden. Das Dataset ist jedoch nicht korrekt an den Grenzen des beabsichtigten Bereichs ausgerichtet. Siehe Bild unten.

Gibt es eine akzeptierte Geotransformation für polare Stereografien, die diese mangelnde Ausrichtung erklären könnte?

Gibt es zum Beispiel einen bestimmten Umrechnungsfaktor oder eine Rotation, die ich verwenden sollte?

Eine Beispieldatei aus dem Datensatz ist hier: "t201113.grd"

Hier ist der Code, den ich derzeit in GDAL verwende

ds = gdal.Open("t201113.grd")
array = ds.ReadAsArray()

x_rotation = 0
y_rotation = 0
xres = 1
yres = -1

llx = -129.8530
lly = 40.0451
ulx = -175.144
uly = 71.385

input_osr = osr.SpatialReference()
input_osr.ImportFromWkt(ds.GetProjection())

wgs84_osr = osr.SpatialReference()
wgs84_osr.ImportFromEPSG(4326)

wgs_to_nps_trans = osr.CoordinateTransformation(wgs84_osr, input_osr)
x, y, z = wgs_to_nps_trans.TransformPoint(llx,lly)

geo_transform = [ x, xres, x_rotation, y, y_rotation, yres ]

ncol = ds.RasterXSize
nrow = ds.RasterYSize

out_driver = gdal.GetDriverByName("HFA")
out_ds = out_driver.Create(t201113.img", ncol, nrow, 1, gdal.GDT_Float32)

out_ds.SetGeoTransform(geo_transform)

out_prj = 'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-110.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",60.0],UNIT["50_Kilometers",50000.0]]'

out_ds.SetProjection(out_prj)

out_ds.GetRasterBand(1).WriteArray(array)
out_ds.GetRasterBand(1).SetNoDataValue(1.70141e+038)
out_ds.FlushCache()
out_ds = None
`

Hier ist auch die Projektionsinformation, wie sie durch die Eingabe definiert ist, dh von "GetProjection ()":

'PROJCS ["North_Pole_Stereographic", GEOGCS ["GCS_WGS_1984", DATUM ["WGS_1984", SPHEROID ["WGS_1984", 6378137.0,298.257223563]], PRIMEM ["Greenwich", 0.0], UNIT ["Degree9.1994. PROJEKTION ["Stereographic"], PARAMETER ["False_Easting", 0.0], PARAMETER ["False_Northing", 0.0], PARAMETER ["Central_Meridian", 0.0], PARAMETER ["Scale_Factor", 1.0], PARAMETER ["Latitude_Of_Origin", 90 ], UNIT ["50_Kilometer", 50000.0]] '

Und die Eingabe GeoTransform:

(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)

Lat, Longs der Gitterkoordinaten werden ebenfalls bereitgestellt, und wenn die Ansicht im projizierten Koordinatensystem wie folgt aussieht. Wenn die Geotransformation durch Koordinaten der unteren linken (gelb) oder oberen rechten (rosa) Koordinate definiert ist, kann ich die Ausdehnung effektiv festlegen, es bleiben jedoch Ausrichtungsprobleme im gesamten Raster.

Bildbeschreibung hier eingeben

jsnider
quelle
Wenn Sie ArcGIS verwenden, wechseln Sie zum stereografischen Nordpol, und legen Sie den Standard parallel auf 60,0 fest. Die stereografische ArcGIS-Implementierung verwendet einen Skalierungsfaktor anstelle einer Standardparallele, da das Projekt an einer beliebigen Stelle zentriert werden kann.
Mkennedy
Thanks @mkennedy - meinst du das Projekt "North Pole Stereographic" (WKID 102018)? Ich habe den Breitengrad des Ursprungs und die Mittelmeridianwerte mithilfe dieser Projektion festgelegt und habe immer noch das gleiche Problem. Vielleicht beziehen Sie sich auf eine andere Projektion?
jsnider
Nein, Sie benötigen eine, bei der die Projektion (Methode) Stereographic_North_Pole ist. Ich glaube nicht, dass wir die genauen PCS haben. Versuchen Sie es mit einer Änderung von 3995 oder 3413.
mkennedy
1
Die Metadaten vermerken: "Die Datei grid_pnt_lls.txt listet die Lat / Longs für jedes x / y auf (0,0 = SW-Ecke des Gitters)." Mit dieser Datei können Sie dieses Raster in ein beliebiges Koordinatensystem umprojizieren und von dort aus fortfahren.
whuber
1
Wo können wir die Vektorebene zum Testen herunterladen?
Farid Cheraghi

Antworten:

2

Zu lange für einen Kommentar, dies ist die Antwort von @ Matej zu begleiten .

  1. Fügen Sie die ".grd" -Daten in ArcGIS hinzu.
  2. Verwenden Sie die Funktion „Raster in anderes Format“ und konvertieren Sie Ihre .grd-Datei in ein ESRII-GRID-Format. Dies ist wichtig, da auf die meisten Rasterfunktionen in ArcGIS nur für dieses Format zugegriffen werden kann. Dies ist entweder der Fall oder es wird normalerweise zu langsam, wenn Sie es in anderen Formaten verwenden.

  3. Da es bereits die Projektionsdatei mit der Datei verknüpft hat. Definieren Sie die Projektion der neuen konvertierten Daten, anstatt sie zu projizieren. ArcToolbox> Data Management Tools> Projektionen und Transformationen> Projektion definieren. Sie können zur vordefinierten polaren Stereo-Grafik-ESRII-Projektion navigieren und prüfen, ob deren Parameter mit den in den Metadaten angegebenen übereinstimmen (nicht). Sie können sie also gemäß @Matej ändern. Nur hier - anstatt zu ändern, erstellen Sie eine neue, die auf der NPS-Projektion mit dem geänderten Mittelmeridian und dem geänderten Breitengrad des Ursprungs basiert, und speichern Sie sie auf der Festplatte. Navigieren Sie dann zur neuen Projektion und verwenden Sie diese, wenn Sie Ihre Projektion definieren. Dies liegt daran, dass Ihre On-the-Fly-Änderung später nicht verfügbar ist, wenn Sie sie zum Festlegen des Koordinatensystems für Ihren Datenrahmen verwenden möchten.

yanes
quelle
1

Ich glaube nicht, dass Sie das Bild neu projizieren müssen. Machen Sie nur Folgendes:

  1. Definieren (Ändern) Sie die Projektion der Grundkarte
  2. Georeferenzieren (verschieben) Sie das Bild an den angegebenen Ort

Beachten Sie, dass sich das Bild (grd) bereits in der StereoGraphic-Projektion für den Nordpol befindet, die nur Aufschluss darüber gibt, wie die Grundkarte angepasst werden soll, die am Bild ausgerichtet wird.

Schritt 1 :

Ändern Sie die ursprüngliche stereografische Projektion des Nordpols (WKID: 102018), um die Breite des Ursprungs und den Mittelmeridian anzupassen:

Bildbeschreibung hier eingeben

Schritt 2:

Referenzieren Sie die grd-Datei, indem Sie die linke untere Ecke auf die angegebene Koordinate (lat, lon) setzen. Wenn Sie die Georeferenzierung aktualisieren, wird die .gdwx-Datei im selben Ordner erstellt. Wenn Sie die SW-Ecke (40.0451, -129.853) zuweisen, sieht der Inhalt der Datei folgendermaßen aus:

50000
0
0
-50000
-1730620.4315
2744092.9724

Bearbeiten: Die obige Weltdatei wurde basierend auf der Zellengröße und der angegebenen Position der SW-Ecke manuell geändert. Die fünfte und sechste Zeile geben die berechnete Position des linken oberen Pixels des Bildes an. Die Position des Bildes hat sich leicht geändert.

Mit den obigen Werten wird das Bild an den angegebenen Ort verschoben und der Maßstab festgelegt.

Und das ist die Ausgabe: Bildbeschreibung hier eingeben

Wenn dies nicht richtig ausgerichtet zu sein scheint, würde ich die angegebenen Koordinaten für die SW-Ecke des Bildes in Frage stellen. Falls Sie Zugriff auf die Koordinaten der NE-Ecke des Bildes haben, können Sie die Transformationsparameter neu berechnen, mit denen das Bild zwischen zwei (oder mehr) Punkten skaliert und gedreht wird.

Matej
quelle
Ist die .gdwx Datei eine Welt - Datei ? In diesem Fall wird in dem verlinkten Wiki-Artikel angegeben, dass die Zeilen 5 und 6 auf das obere linke Pixel verweisen . Schlagen Sie vor, es auf den Südwestpixel (unten links) einzustellen?
Kirk Kuykendall
Nein. Ich habe nur die Position der SW-Ecke wie in der Readme-Datei angegeben angegeben. Die Struktur der Datei sieht aus wie die Weltdatei. Es wurde mit dem Georeferenzierungs-Tool in ArcMap generiert, das möglicherweise die Position der NW-Ecke berechnet und gespeichert hat - noch nicht überprüft.
Matej
1
Ja, ich habe es jetzt überprüft. Der gespeicherte Ort in der .gdwx ist die obere linke Ecke.
Matej,