Геотрансформация для полярной стереографии?

16

В настоящее время я работаю над тем, чтобы импортировать климатические данные CANGRID (предоставленные как Surfer Grid ascii, файлы .grd) в ArcGIS. Сетка размером 95 строк на 125 столбцов. Метаданные обеспечивают широту / долготу происхождения (нижний левый угол), размер ячейки (50 км), а также проекцию нот в виде полярной стереографии с центральным меридианом (110 градусов западной долготы) и широтой происхождения (60 градусов северной широты).

После первой попытки безуспешно преобразовать .grd в растры как .ascii и .flt, мне удалось использовать GDAL для установки экстента и проекции, однако набор данных неправильно совмещается с границами предполагаемой области. Смотрите изображение ниже.

Существует ли общепринятая геотрансформация для полярной стереографии, которая могла бы объяснить это отсутствие выравнивания?

Например, есть ли конкретный коэффициент преобразования или ротация, который я должен использовать?

Пример файла из набора данных находится здесь: "t201113.grd"

Вот код, который я сейчас использую в GDAL

ds = gdal.Open("t201113.grd")
array = ds.ReadAsArray()

x_rotation = 0
y_rotation = 0
xres = 1
yres = -1

llx = -129.8530
lly = 40.0451
ulx = -175.144
uly = 71.385

input_osr = osr.SpatialReference()
input_osr.ImportFromWkt(ds.GetProjection())

wgs84_osr = osr.SpatialReference()
wgs84_osr.ImportFromEPSG(4326)

wgs_to_nps_trans = osr.CoordinateTransformation(wgs84_osr, input_osr)
x, y, z = wgs_to_nps_trans.TransformPoint(llx,lly)

geo_transform = [ x, xres, x_rotation, y, y_rotation, yres ]

ncol = ds.RasterXSize
nrow = ds.RasterYSize

out_driver = gdal.GetDriverByName("HFA")
out_ds = out_driver.Create(t201113.img", ncol, nrow, 1, gdal.GDT_Float32)

out_ds.SetGeoTransform(geo_transform)

out_prj = 'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-110.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",60.0],UNIT["50_Kilometers",50000.0]]'

out_ds.SetProjection(out_prj)

out_ds.GetRasterBand(1).WriteArray(array)
out_ds.GetRasterBand(1).SetNoDataValue(1.70141e+038)
out_ds.FlushCache()
out_ds = None
`

Кроме того, вот информация о проекции, как определено вводом, то есть из «GetProjection ()»:

«PROJCS [ "North_Pole_Stereographic", GEOGCS [ "GCS_WGS_1984", НУЛЕВОЙ [ "WGS_1984", сфероида [ "WGS_1984", 6378137.0,298.257223563]], PRIMEM [ "Гринвич", 0,0], блок [ "Степень", 0,0174532925199433]], ПРОЕЦИРОВАНИЯ [ "стереографическая"], параметр [ "False_Easting", 0,0], параметр [ "False_Northing", 0,0], параметр [ "Central_Meridian", 0,0], параметр [ "Scale_Factor", 1,0], параметр [ "Latitude_Of_Origin", 90.0 ], [БЛОК "50_Kilometers", 50000,0]]»

И вход GeoTransform:

(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)

Широта, длинные координаты сетки также предоставляются, и при просмотре в проецируемой системе координат выглядят как показано ниже. Когда геотрансформация определяется координатами нижнего левого (желтого) или верхнего правого (розового) кордината, я могу эффективно установить экстент, но остаются проблемы с выравниванием по всему растру.

введите описание изображения здесь

jsnider
источник
Если вы используете ArcGIS, переключитесь на Стереографический Северный полюс и установите стандартную параллель на 60.0. Стереографическая реализация ArcGIS использует масштабный коэффициент, а не стандартную параллель, потому что проект может быть центрирован где угодно.
Mkennedy
Спасибо @mkennedy - Вы имеете в виду проект "Стереографический Северный полюс" (WKID 102018)? Я установил широту происхождения и значения центрального меридиана, используя эту проекцию, и все еще имею ту же проблему. Возможно, вы имеете в виду другую проекцию?
jsnider
Нет, вам нужен такой, где проекция (метод) - Stereographic_North_Pole. Я не думаю, что у нас есть точный PCS; попробуйте изменить с 3995 или 3413.
Mkennedy
1
В метаданных отмечается, что «В файле grid_pnt_lls.txt перечислены широты / долготы для каждого x / y (0,0 = угол сетки SW)». Имея этот файл в руках, вы можете перепроецировать эту сетку в любую нужную вам систему координат и идти дальше.
whuber
1
Где можно скачать векторный слой для тестирования?
Фарид Чераги

Ответы:

2

Слишком долго для комментария, это должно сопровождать ответ @ Matej .

  1. Добавьте данные «.grd» в ArcGIS.
  2. Используйте функцию «Растр в другой формат» и преобразуйте ваш файл .grd в формат ESRII GRID. Это важно, потому что большинство растровых функций в ArcGIS доступны только для этого формата, либо он обычно слишком медленный, когда вы используете его в других форматах.

  3. Так как он уже имеет проекционный файл, связанный с файлом. Вместо того, чтобы проецировать новые преобразованные данные, определите их проекцию. ArcToolbox> Инструменты управления данными> Проекции и преобразования> Определить проекцию. Вы можете перейти к предварительно определенной полярно-стереографической графической проекции ESRII и посмотреть, соответствуют ли ее параметры параметрам, указанным в метаданных (это не так), поэтому вы можете изменить его согласно @Matej. Только здесь - вместо изменения создайте новый на основе проекции NPS с измененным центральным меридианом и широтой координат и сохраните его на диске, затем перейдите к новой проекции и используйте ее при определении проекции. Это потому, что ваша модификация на лету не будет доступна позже, когда вы захотите использовать ее для установки системы координат для вашего фрейма данных,

Yanes
источник
1

Я не думаю, что вам нужно перепроектировать изображение. Только сделайте следующее:

  1. определить (изменить) проекцию базовой карты
  2. привязка (смещение) изображения к указанному месту

Обратите внимание, что изображение (grd) уже находится в стерео-графической проекции Северного полюса, которая дает вам только представление о том, как настроить базовую карту, которая будет выровнена с изображением.

Шаг 1 :

Измените исходную стереографическую проекцию Северного полюса (WKID: 102018), чтобы настроить широту происхождения и центральный меридиан:

введите описание изображения здесь

Шаг 2:

При помощи геопривязки к файлу grd установите нижний левый угол для указанной координаты (широта, долгота). При обновлении географической привязки файл .gdwx создается в той же папке. При назначении угла SW для (40.0451, -129.853) содержимое файла выглядит следующим образом:

50000
0
0
-50000
-1730620.4315
2744092.9724

редактировать: файл мира выше был изменен вручную на основе размера ячейки и предоставленного местоположения угла SW - 5-я и 6-я строки представляют расчетное местоположение верхнего левого пикселя изображения. Положение изображения немного изменилось.

Вышеуказанные значения помещают (смещают) изображение в указанное место и определяют масштаб.

И это вывод: введите описание изображения здесь

Если это не похоже на правильное выравнивание, я бы поставил под сомнение предоставленные координаты для угла SW изображения. Если у вас есть доступ к координатам, т. Е. К углу NE изображения, вы можете пересчитать параметры преобразования, которые бы масштабировали и поворачивали изображение между двумя (или более) точками.

Матей
источник
Является ли файл .gdwx файлом мира ? Если это так, то в связанной вики статье говорится, что строки 5 и 6 ссылаются на верхний левый пиксель. Вы предлагаете установить его на юго-западный (нижний левый) пиксель?
Кирк Куйкендалл
Нет. Я указал только местоположение угла SW, как отмечено в файле readme. Структура файла выглядит как мировой файл. Он был сгенерирован с использованием инструмента Georeferencing в ArcMap, который, возможно, рассчитал и сохранил местоположение угла NW - еще не проверил.
Матей
1
Да, я проверил это сейчас. Сохраненное местоположение в .gdwx - верхний левый угол.
Матей