Mais

Usando Proj.4 / GDAL para converter lat, lon em i, j

Usando Proj.4 / GDAL para converter lat, lon em i, j


Estou trabalhando com um conjunto de dados no formato GRIB2 e preciso determinar o índice de células da grade (eu j) correspondendo a um determinado lat, lon. Tenho a definição de projeção na forma de projparams da seguinte forma:

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

o que me dá os seguintes srs:+ unidades = m + a = 6371229 + b = 6371229 + lon_0 = 265,0 + proj = lcc + lat_2 = 25,0 + lat_1 = 25,0 + lat_0 = 25,0

E eu sei dx, dy, número de pontos de grade em cada direção etc.

Então, eu estava pensando em usar Proj.4 (porque tem vínculos JavaScript) (e / ou GDAL) para convertê-lo em i, j. Acho que li em algum lugar que se poderia fazer dessa forma, definindo algum tipo de projeção de mapa métrico equidistante que usariaeu jComox, ycoordena e faz uma transformação de coordenada simples, mas estou lutando para definir como realmente essa projeção em termos de proj.4.


Aqui está como o utilitário wgrib2 faz isso. O algoritmo:

  1. Inicialize 2 projeções: um x, y chamadopj_gride um lat, lon chamadopj_latlon
  2. Dado lat, lon da primeira célula da grade, encontre suas coordenadas x, y (x_0, y_0)
  3. Dado lat, lon de um ponto de interesse, encontre suas coordenadas x, y (x, y)
  4. Dados os espaçamentos da grade dx, dy, encontre os índices da grade do ponto de interesseeu jComo:

    i = (x-x_0) / dx + 1

    j = (y-y_0) / dy + 1

Nota: os índices são baseados em 1. Abaixo está um exemplo de implementação em Python usando pyproj. A grade usada é a grade NAM218

import pyproj # no meu caso projparams é obtido usando 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 # + unidades = 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 # espaçamento da grade em m dx = 12191 dy = 12191 # primeiro ponto da grade, lat / lon em graus lat1 = 12,19 lon1 = 226,541 x_0, y_0 = pyproj.transform (pj_latlon, pj_grid, lon1, lat1, radianos = Falso) imprimir x_0, y_0 # -4226106.99692 -832698.261018 lat = 51 lontransform (y = pyproj. pj_latlon, pj_grid, lon, lat, radianos = False) # estes são índices baseados em 1 i = (x - x_0) / dx + 1 j = (y - y_0) / dy + 1

Assista o vídeo: QGIS Basic #102: Lat Lon Tools Part 1