Как проецировать координаты x, y в широты / долготы в файле netcdf - PullRequest
1 голос
/ 02 июля 2019

Я скачал поле скоростей ледового щита Гренландии с веб-сайта CCI в виде файла NetCDF.Однако проекция задается как (см. Ниже, где x находится в диапазоне от [-639750,855750] до y [-655750, -3355750])

Как проецировать эти данные в фактические координаты широты / долготы вфайл NetCDF?Уже спасибо!Для интересующихся: файл можно скачать здесь: http://products.esa -icesheets-cci.org / products / downloadlist / IV /

Variables:
    crs                                
           Size:       1x1
           Dimensions: 
           Datatype:   int32
           Attributes:
                       grid_mapping_name                     = 'polar_stereographic'
                       standard_parallel                     = 70
                       straight_vertical_longitude_from_pole = -45
                       false_easting                         = 0
                       false_northing                        = 0
                       unit                                  = 'meter'
                       latitude_of_projection_origin         = 90
                       spatial_ref                           = 'PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3413"]]'
    y                                  
           Size:       5401x1
           Dimensions: y
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'Y'
                       long_name     = 'y coordinate of projection'
                       standard_name = 'projection_y_coordinate'
    x                                  
           Size:       2992x1
           Dimensions: x
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'X'
                       long_name     = 'x coordinate of projection'
                       standard_name = 'projection_x_coordinate'

1 Ответ

2 голосов
/ 05 июля 2019

Если вы хотите преобразовать всю сетку из ее собственных полярных стереографических координат в географическую (долгота по широте) сетку, вы, вероятно, захотите использовать такой инструмент, как gdalwarp . Я не думаю, что вы задаете этот вопрос.

Если я правильно читаю ваш вопрос, вы хотите выбрать точки из файла и найти их как пары координат lon / lat. Я предполагаю, что вы знаете, как получить местоположение в виде пары XY из вашего файла netCDF вместе со значениями скорости в этом месте. Я также предполагаю, что вы делаете это в Python, так как вы поставили этот тег на этот вопрос.

Когда у вас есть пара XY, вам просто нужна функция (с кучей параметров), чтобы преобразовать ее в lon / lat. Вы можете найти эту функцию в модуле pyproj .

Pyproj оборачивает библиотеку proj4 C, которая очень широко используется для преобразований системы координат. Если у вас есть пара XY в проецируемых координатах и ​​вы знаете определение проецируемой системы координат, вы можете использовать функцию преобразования pyproj, например:

import pyproj

# Output coordinates are in WGS 84 longitude and latitude
projOut = pyproj.Proj(init='epsg:4326')

# Input coordinates are in meters on the Polar Stereographic 
# projection given in the netCDF file
projIn = pyproj.Proj('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 
    +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ',
    preserve_units=True)

# here is a coordinate pair near the middle of your data set
x, y = 0.0, -2000000

# transform x,y to lon/lat
lon, lat = pyproj.transform(projIn, projOut, x, y)

# answer: lon = -45.0; lat = 71.6886

... и вот, пожалуйста. Обратите внимание, что выходная долгота равна -45.0, что должно дать вам приятное ощущение тепла, поскольку координата X входного сигнала была 0, а -45.0 - центральный меридиан проекции набора данных. Если вы хотите получить ответ в радианах, а не в градусах, установите radians kwarg в функции преобразования на True.

Теперь самое сложное, что вы и делаете в первую очередь - определяете projIn и projOut, которые используются в качестве аргументов для функции преобразования. Они находятся во входных и выходных системах координат для преобразования. Это объекты Proj, и они содержат множество параметров для уравнений преобразования системы координат. Разработчики proj4 заключили их в аккуратный набор функций, а разработчики pyproj обернули вокруг них красивую оболочку python, так что вам и мне не нужно отслеживать все детали. Я буду благодарен им за все оставшиеся дни.

Выходная система координат тривиальна

projOut = pyproj.Proj(init='epsg:4326')

Библиотека pyproj может создавать объект Proj из кода EPSG. 4326 - это код EPSG для WGS 84 lon / lat. Готово.

Установка projIn сложнее, потому что ваш файл netCDF определяет свою систему координат с помощью строки WKT, которую (я почти уверен) не может быть прочитан напрямую proj4 или pyproj. Однако pyproj.Proj () будет принимать строку параметра proj4 в качестве аргумента. Я уже дал вам тот, который вам нужен для этой операции, так что вы можете просто взять мой за это, что это

+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

является эквивалентом этого (который копируется непосредственно из вашего файла netCDF):

PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
  GEOGCS["WGS 84",
    DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Polar_Stereographic"],
  PARAMETER["latitude_of_origin",70],
  PARAMETER["central_meridian",-45],
  PARAMETER["scale_factor",1],
  PARAMETER["false_easting",0],
  PARAMETER["false_northing",0],
  UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
AUTHORITY["EPSG","3413"]]'

Если вы хотите сделать это в более общем плане, вам понадобится другой модуль для преобразования определений системы координат WKT в строки параметров proj4. Одним из таких модулей является osgeo.osr, и в этом блоге есть пример программы, в которой показано, как конвертировать .

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...