Geotransformasi untuk stereografi kutub?

16

Saat ini saya sedang bekerja untuk mengimpor data iklim CANGRID (disediakan sebagai Surfer Grid ascii, file ".grd") ke ArcGIS. Grid berukuran 95 baris dengan 125 kolom. Metadata menyediakan lat / lon asal (sudut kiri bawah), ukuran sel (50 km) serta catatan proyeksi sebagai stereografi polar dengan meridian pusat (110 derajat W) dan garis lintang asal (60 derajat N).

Setelah pertama kali mencoba mengubah .grd menjadi raster sebagai .ascii dan .flt tidak berhasil, saya berhasil menggunakan GDAL untuk menetapkan tingkat dan proyeksi, namun set data tidak benar sesuai dengan batas-batas area yang dimaksud. Lihat gambar di bawah.

Apakah ada geotransformasi yang diterima untuk stereografi polar yang dapat menjelaskan kurangnya penyelarasan ini?

Misalnya, apakah ada faktor konversi atau rotasi tertentu yang harus saya gunakan?

Contoh file dari dataset ada di sini: "t201113.grd"

Berikut kode yang saya gunakan saat ini di 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
`

Juga, inilah info proyeksi, seperti yang didefinisikan oleh input, yaitu dari "GetProjection ()":

'PROJCS ["North_Pole_Stereographic", GEOGCS ["GCS_WGS_1984", DATUM ["WGS_1984", SPHEROID ["WGS_1984", 6378137.0.2998.257223563]], PRIMEM ["Greenwich", 0,0], 0,01253 "0,334,324" PROYEKSI ["Stereografis"], PARAMETER ["False_Easting", 0,0], PARAMETER ["False_Northing", 0,0], PARAMETER ["Central_Meridian", 0,0], PARAMETER ["Scale_Factor", 1.0], PARAMETER ["Latitude_Of_Origin, ".0.0, 90.0, 90.0) ], UNIT ["50_Kilometer", 50000.0]] '

Dan input GeoTransform:

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

Lat, panjang koordinat grid juga disediakan, dan ketika tampilan dalam sistem koordinat yang diproyeksikan terlihat seperti di bawah ini. Ketika geotransform didefinisikan oleh koordinat kordinat kiri bawah (kuning) atau kanan atas (merah muda), saya dapat secara efektif mengatur luasnya, tetapi masih ada masalah penyelarasan di seluruh raster.

masukkan deskripsi gambar di sini

jsnider
sumber
Jika Anda menggunakan ArcGIS, beralihlah ke Kutub Utara Stereografis dan atur paralel standar ke 60.0. Implementasi stereografik ArcGIS menggunakan faktor skala daripada standar paralel karena proj dapat dipusatkan di mana saja.
mkennedy
Terima kasih @kenkeny - maksud Anda proj "Stereografi Kutub Utara" (WKID 102018)? Saya telah menetapkan garis lintang asal dan nilai meridian pusat menggunakan proyeksi ini dan masih memiliki masalah yang sama. Mungkin Anda mengacu pada proyeksi lain?
jsnider
Tidak, Anda memerlukannya di mana proyeksi (metode) adalah Stereographic_North_Pole. Saya tidak berpikir kita memiliki PCS yang tepat; coba modifikasi dari 3995 atau 3413.
mkennedy
1
Metadata mencatat bahwa "File grid_pnt_lls.txt mencantumkan lat / long untuk setiap x / y (0,0 = sudut SW dari kisi-kisi)." Dengan file ini di tangan Anda dapat memproyeksi ulang grid ini ke sistem koordinat yang Anda inginkan dan melanjutkan dari sana.
whuber
1
Di mana kita bisa mengunduh layer vektor untuk diuji?
Farid Cheraghi

Jawaban:

2

Terlalu lama untuk berkomentar, ini untuk menemani jawaban @ Matej .

  1. Tambahkan data ".grd" ke dalam ArcGIS.
  2. Gunakan fungsi “Raster to Other Format” dan konversikan file .grd Anda menjadi format ESRII GRID. Ini penting karena sebagian besar fungsi raster di ArcGIS hanya dapat diakses untuk format ini, baik itu atau biasanya terlalu lambat ketika Anda menggunakannya pada format lain.

  3. Karena sudah memiliki file proyeksi yang terkait dengan file tersebut. Daripada memproyeksikan data yang dikonversi baru, tentukan proyeksinya. ArcToolbox> Alat Manajemen Data> Proyeksi dan Transformasi> Tetapkan Proyeksi. Anda dapat menavigasi ke proyeksi ESRII grafik stereo polar yang telah ditentukan sebelumnya dan melihat apakah parameternya cocok dengan yang diberikan dalam metadata (tidak), sehingga Anda dapat memodifikasinya sesuai @Matej. Hanya di sini - alih-alih memodifikasi, buat yang baru berdasarkan proyeksi NPS dengan meridian pusat dan Latitude asal diubah dan simpan di disk, lalu navigasikan ke proyeksi baru dan gunakan itu saat menentukan proyeksi Anda. Ini karena modifikasi on the fly Anda tidak akan tersedia nanti ketika Anda ingin menggunakannya untuk mengatur sistem koordinat untuk bingkai data Anda,

yanes
sumber
1

Saya tidak berpikir Anda perlu memproyeksi ulang gambar. Hanya lakukan hal berikut:

  1. tentukan (modifikasi) proyeksi basemap
  2. georeferensi (menggeser) gambar ke lokasi yang ditentukan

Perhatikan bahwa gambar (grd) sudah ada dalam proyeksi StereoGraphic Kutub Utara, yang hanya memberi Anda indikasi cara menyesuaikan basemap yang akan disejajarkan dengan gambar.

Langkah 1 :

Ubah proyeksi Stereografis Kutub Utara asli (WKID: 102018) untuk menyesuaikan Latitude of Origin dan Central Meridian:

masukkan deskripsi gambar di sini

Langkah 2:

Georeferensi file grd dengan mengatur sudut kiri bawah ke koordinat yang ditentukan (lat, lon). Ketika Anda memperbarui georeferensi, file .gdwx dibuat di folder yang sama. Saat menetapkan sudut SW ke (40.0451, -129.853) konten file terlihat sebagai berikut:

50000
0
0
-50000
-1730620.4315
2744092.9724

sunting: file dunia di atas telah dimodifikasi secara manual berdasarkan ukuran sel dan lokasi sudut SW yang disediakan - garis ke-5 dan ke-6 merupakan lokasi terhitung dari piksel kiri atas gambar. Posisi gambar sedikit berubah.

Nilai-nilai di atas menempatkan (menggeser) gambar ke lokasi yang ditentukan dan menentukan skala.

Dan inilah hasilnya: masukkan deskripsi gambar di sini

Jika ini tampaknya tidak selaras dengan benar, saya akan mempertanyakan koordinat yang disediakan untuk sudut SW gambar. Jika Anda memiliki akses ke koordinat yaitu sudut NE gambar, Anda bisa menghitung ulang parameter transformasi yang akan skala dan memutar gambar antara dua (atau lebih) titik.

Matej
sumber
Apakah file .gdwx adalah file dunia ? Jika demikian, maka artikel wiki yang tertaut mengatakan baris 5 & 6 merujuk piksel kiri atas . Apakah Anda menyarankan pengaturan ke piksel barat daya (kiri bawah)?
Kirk Kuykendall
Tidak. Saya hanya menentukan lokasi sudut SW seperti yang tercantum dalam file readme. Struktur file terlihat seperti file dunia. Itu dihasilkan menggunakan alat Georeferencing di ArcMap yang mungkin telah menghitung dan menyimpan lokasi sudut NW - belum diperiksa.
Matej
1
Ya, saya memeriksanya sekarang. Lokasi yang disimpan di .gdwx adalah sudut kiri atas.
Matej