Suite

Utilisation de Proj.4/GDAL pour convertir lat,lon en i,j

Utilisation de Proj.4/GDAL pour convertir lat,lon en i,j


Je travaille avec un ensemble de données au format GRIB2 et j'ai besoin de déterminer l'index des cellules de la grille (je, j) correspondant à une lat,lon donnée. J'ai la définition de la projection sous la forme de projparams comme suit :

{'a' : 6371229, 'b' : 6371229, 'lat_0' : 25,0, 'lat_1' : 25,0, 'lat_2' : 25,0, 'lon_0' : 265,0, 'proj' : 'lcc'}

ce qui me donne les srs suivants :+units=m +a=6371229 +b=6371229 +lon_0=265.0 +proj=lcc +lat_2=25.0 +lat_1=25.0 +lat_0=25.0

Et je connais dx, dy, le nombre de points de grille dans chaque direction, etc.

Je pensais donc utiliser Proj.4 (car il a des liaisons JavaScript) (et/ou GDAL) pour le convertir en i,j. Je pense avoir lu quelque part que l'on pourrait procéder de cette façon en définissant une sorte de projection cartographique métrique équidistante qui utiliseraitje, jcommex, ycoordonnées et faire une simple transformation de coordonnées, mais j'ai du mal à définir réellement une telle projection en termes de proj.4.


Voici comment l'utilitaire wgrib2 le fait. L'algorithme :

  1. Initialiser 2 projections : une x,y appeléepj_gridet un lat,lon appelépj_latlon
  2. Étant donné lat,lon de la première cellule de la grille, trouvez ses coordonnées x,y (x_0, y_0)
  3. Étant donné lat,lon d'un point d'intérêt, trouver ses coordonnées x,y (x, y)
  4. Étant donné les espacements de grille dx, dy, trouver les indices de grille du point d'intérêtje, jcomme:

    i = (x-x_0)/dx + 1

    j = (y-y_0)/j + 1

Remarque : les indices sont basés sur 1. Vous trouverez ci-dessous un exemple d'implémentation en Python à l'aide de pyproj. La grille utilisée est la grille NAM218

import pyproj # dans mon cas, projparams est obtenu en utilisant pygrib 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 # espacement de la grille en m dx = 12191 dy = 12191 # premier point de grille, lat/lon en degrés lat1 = 12.19 lon1 = 226.541 x_0, y_0 = pyproj.transform(pj_latlon, pj_grid, lon1, lat1, radians=False) print x_0, y_0 # -4226106.99692 -832698.261018 lat = 51 lon = -114 x, y = (pyproj. pj_latlon, pj_grid, lon, lat, radians=False) # ce sont des indices basés sur 1 i = (x - x_0) / dx + 1 j = (y - y_0) / dy + 1

Voir la vidéo: How to find Latitude Longitude from Topographic Map