Если вы хотите преобразовать всю сетку из ее собственных полярных стереографических координат в географическую (долгота по широте) сетку, вы, вероятно, захотите использовать такой инструмент, как 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, и в этом блоге есть пример программы, в которой показано, как конвертировать .