Jak przekonwertować współrzędne afiniczne na lat / lng?

14

Jestem wyjątkowo nowy w GIS.

Używam gdaldo czytania na mapie użytkowania terenu / pokrycia terenu i muszę wybrać szerokość / długość niektórych rodzajów pokrycia terenu, aby zindeksować do innego zestawu danych, który jest wyrażony tylko w szerokości / długości. Niestety nie rozumiem formy współrzędnych xiy podanych mi z geotransformy, a konkretnie tych originXi originYponiżej:

geotransform = dataset.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]

Drukowanie tych wartości daje mi współrzędne jak (447466.693808, 4952570.40529). Jak odnoszą się one do pierwotnej szerokości i długości geograficznej?

Edytować:

Oto prosty przykład python, który dał mi to, czego szukałem:

srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)

Skradzione z: tolatlong.py

gdal coordinate-system

— Bogaty
źródło

Wygląda na to, że Twoje dane są wyświetlane (np. UTM) i musisz wiedzieć, czym jest projekcja, aby „odrzucić” ją z powrotem do współrzędnych długich / długich.

@Dan Dzięki, więc wiem, że mogę uzyskać projekcję za pomocą a dataset.GetProjectionRef()i dowiedzieć się, że używam „UTM Zone 10”, ale co wtedy? Poszukuję metod takich jak „nieprojekt”, ale zbliżam się do zera.

— Bogaty

przepraszam za użycie terminu „nieprojekt” (w cudzysłowach), ponieważ dane w stopniach dziesiętnych nie są rzutowane, ale jeśli chcesz odzyskać rzutowane dane z powrotem do stopni dziesiętnych z dowolnej projekcji, musisz (zwróć uwagę na cytaty) „projekt” wraca do układu współrzędnych geograficznych, czyli danych w stopniach dziesiętnych.

2

Ten (nowszy) wątek zawiera wyraźny przykład i inne rozwiązanie: gis.stackexchange.com/questions/8430/…

— whuber

10

gdal_translate przerzuci twoje dane z dowolnej projekcji na dowolną inną (w tym przypadku chcesz EPSG: 4326), używając:

gdal_translate -a_srs epsg:4326 srcfile new_file

lub możesz użyć gdaltrasform do konwersji punktów (i jestem pewien, że możesz uzyskać do tego dostęp również z Python (?))

— Ian Turton
źródło

(+1) Nie wiem, dlaczego to zostało przegłosowane, Ian, ponieważ wydaje mi się, że to jakaś operacja, która będzie potrzebna po zastosowaniu transformacji afinicznej.

— whuber

gdal_translate nie będzie działał tak, jak go masz (a raczej po prostu przypisałby EPSG: 4326 do danych), musisz wypaczać dane za pomocą czegoś takiego: gdalwarp -s_srs EPSG: 32610 -t_srs EPSG: 4326 src dest Jeśli dane mają już metadane proj, możesz porzucić parametr -s_srs.

— MerseyViking,

Może właśnie o to mi chodzi. Czy „epsg: 4326” jest przekształceniem w globalny język lat / lng? Wydaje mi się, że widzę niektóre klasy, w których mogę przedstawić projekcję, myślałem „WGS: 84”, ale nie znam różnicy.

— Bogaty

1

@ Rich Kliknij znacznik „układ współrzędnych” w swoim pytaniu, posortuj wynikową stronę według głosów i zacznij czytać. Odpowiedzi na wiele pytań udzielisz w ciągu kilku minut.

— whuber

7

Geotransform jest udokumentowany na https://gdal.org/user/raster_data_model.html . Chodzi o to, że bierzesz współrzędne (x, y) prosto z zestawu danych, stosujesz transformację liniową, aby uzyskać (u, v) z

u = a*x + b*y
v = c*x + d*y

(można przyjąć, że jest to definicja transformacji liniowej), a następnie przesunąć początek, dodając geotransform [0] do u i geotransform [3] do v. To daje „transformację afiniczną” (x, y). Naprawdę ma na celu obracanie, zmianę skali, być może trochę poprawienie niektórych błędów skosu i ponowne ustawienie współrzędnych specyficznych dla danych (x, y), aby pasowały do znanego układu współrzędnych. Wynik ma dawać współrzędne rzutowane. Oznacza to po prostu, że procedura matematyczna przyjmuje (długość, szerokość) i zamienia je na znane współrzędne: nazywa się to „rzutowaniem”. „Unprojecting” robi coś przeciwnego; więc jeśli wiesz, która projekcja jest potrzebna, zastosuj ją do współrzędnych przekształconych współrzędnych (x, y) aby uzyskać szerokość i długość geograficzną.

Nawiasem mówiąc, wartości stałych a, b, c, d są podane przez 1, 2, 4 i 5 wpisów w tablicy geotransformacyjnej.

— Whuber
źródło

1

Przeczytałem sekcję geotransform podanego linku, ale nie sądzę, że o to mi chodzi. Przepraszam, jeśli jestem tępy, ale czy to nie opisuje, jak rzutować indeksy xiy (0 na XSize lub YSize) na rzutowane współrzędne? Jestem ciekawy, jak tłumaczysz rzutowane współrzędne na szerokość i długość geograficzną. Jako prosty przykład geotransforma definiuje początek x / y lewego górnego rogu rastra w rzutowanych współrzędnych (447466.693808, 4952570.40529). Jak zamieniam te współrzędne z powrotem na lat / lng (coś w rodzaju lat: 44, lng: -122).

— Bogaty

@Rich Nie, ten link nie opisuje projekcji. Opisuje jedynie zmianę współrzędnych między dwiema mapami, a nie między światem a mapą. Unprojection obraca afiniczna transformacji współrzędnych do (szer). Zdajesz nie chce napisać kod, że jeśli możesz pomóc! Unprojection prawie zawsze polega na określeniu potrzebnej projekcji i wywołaniu odpowiedniego oprogramowania, które wykona zadanie za Ciebie (jak sugeruje odpowiedź @ iant).

— whuber

Link jest zepsuty. @ czy mam rację, że transformacja afiniczna (z GetGeoTransform GDALa) jest transformacją między pikselami a współrzędnymi CRS? Na przykład, jeśli dane są w projekcji GDA, punkty należy przekształcić z WGS84 / lat / lon do GDAXX przy użyciu np. CoordianteTransform, a następnie w piksele przy użyciu afinicznego GeoTransform? Ten dwuetapowy proces bardzo mnie potknął i podejrzewam, że również powoduje problemy z OP. Nie pomogło, że xarray / rasterio wydaje się mieć błąd, który psuje wartości GeoTransform

— naught101

@naught Znalazłem nowy adres URL i zaktualizowałem swój post.

— whuber

0

Możesz użyć następujących opcji:

import gdal, numpy as n
def coord(file):
    padfTransform = file.GetGeoTransform()
    indices = n.indices(file.ReadAsArray().shape)
    xp = padfTransform[0] + indices[1]*padfTransform[1] + indices[1]*padfTransform[2]   
    yp = padfTransform[3] + indices[0]*padfTransform[4] + indices[0]*padfTransform[5]  
    return xp,yp
file = gdal.Open('Yourfile.tif') # A GeoTiff file
x,y = coord(file)

koordyn zwróci długość (x) i szerokość (y) wszystkich pikseli. Pamiętaj, że współrzędne znajdują się w lewym rogu piksela

— Ishan Tomar
źródło

Czy nie powinno to oznaczać [1] * padfTransform [1] + indeksów [0] * padfTransform [2] i odwrotnie dla YP?

— CMCDragonkai