Mehr

Verwenden von Proj.4/GDAL zum Konvertieren von lat,lon in i,j

Verwenden von Proj.4/GDAL zum Konvertieren von lat,lon in i,j


Ich arbeite mit einem Datensatz im GRIB2-Format und muss den Rasterzellenindex (ich, j) entsprechend einem gegebenen lat,lon. Ich habe die Projektionsdefinition in Form von projparams wie folgt:

{'a': 6371229, 'b': 6371229, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}

was mir folgende srs gibt:+Einheiten=m +a=6371229 +b=6371229 +lon_0=265.0 +proj=lcc +lat_2=25,0 +lat_1=25,0 +lat_0=25,0

Und ich kenne dx, dy, Anzahl der Gitterpunkte in jede Richtung usw.

Also dachte ich daran, Proj.4 (weil es JavaScript-Bindungen hat) (und/oder GDAL) zu verwenden, um es in i,j zu konvertieren. Ich glaube, ich habe irgendwo gelesen, dass man es auf diese Weise machen könnte, indem man eine Art äquidistante metrische Kartenprojektion definiert, die verwenden würdeich, jwiex, yKoordinaten und mache eine einfache Koordinatentransformation, aber ich kämpfe damit, wie man eine solche Projektion tatsächlich in Proj.4-Begriffen definieren kann.


So macht es das Dienstprogramm wgrib2. Der Algorithmus:

  1. Initialisieren Sie 2 Projektionen: eine XY-Projektion namenspj_gridund ein lat, lon nannte manpj_latlon
  2. Gegeben lat,lon der ersten Gitterzelle, bestimme ihre x,y-Koordinaten (x_0, y_0)
  3. Gegeben lat,lon eines Point of Interest, bestimme seine x,y-Koordinaten (x,y)
  4. Bestimmen Sie mit den Gitterabständen dx, dy die Gitterindizes des Point of Interestich, jwie:

    i = (x-x_0)/dx + 1

    j = (y-y_0)/dy + 1

Hinweis: Die Indizes basieren auf 1. Unten sehen Sie eine Beispielimplementierung in Python mit pyproj. Das verwendete Raster ist NAM218 Raster

import pyproj # in meinem Fall wird projparams mit pygrib abgerufen projparams = { 'a': 6371229, 'b': 6371229, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc' } pj_grid = pyproj.Proj(projparams=projparams) print pj_grid.srs # +units=m +a=6371229 +b=6371229 +lon_0=265.0 +proj=lcc +lat_2=25.0 +lat_1= 25.0 +lat_0=25.0 pj_latlon = pj_grid.to_latlong() print pj_latlon.srs # +proj=latlong +a=6371229 +b=6371229 # Rasterabstand in m dx = 12191 dy = 12191 # erster Rasterpunkt, Breite/Länge in Grad lat1 = 12.19 lon1 = 226.541 x_0, y_0 = pyproj.transform(pj_latlon, pj_grid, lon1, lat1, Radiant=False) print x_0, y_0 # -4226106.99692 -832698.261018 lat = 51 lon = -114 x, y = py pj_latlon, pj_grid, lon, lat, radians=False) # dies sind 1-basierte Indizes i = (x - x_0) / dx + 1 j = (y - y_0) / dy + 1

Schau das Video: Convert UTM Coordinates to Latitude Longitude. ENZ to Latitude u0026 Longitude in Global Mapper