Przeprojektowanie WGS 1984 Web Mercator (EPSG: 3857) w Pythonie za pomocą GDAL

17

Ponownie projektuję rastry w Pythonie za pomocą GDAL. Muszę rzutować kilka tiffów ze współrzędnych geograficznych WGS 84 na WGS 1984 Web Mercator (Sfera Pomocnicza), aby użyć ich później w Openlayers razem z OpenStreetMap i być może mapami Google. Korzystam z Pythona 2.7.5 i GDAL 1.10.1 stąd i przekształcam współrzędne za pomocą porad z tego miejsca (mój kod jest poniżej). Krótko mówiąc, zaimportowałem osgeo.osr i użyłem ImportFromEPSG (kod) i CoordinateTransformation (z, do) .

Najpierw próbowałem EPSG (32629), który jest strefą UTM 29 i otrzymałem ten rzutowany raster (mniej więcej dobrze), więc kod wydaje się poprawny: utm Potem użyłem EPSG (3857), ponieważ przeczytałem to i to pytania i znalazłem że jest to poprawny aktualny poprawny kod . Ale raster jest tworzony bez odniesienia przestrzennego. Jest daleko w ramce danych WGS 84 (ale będzie dobrze, jeśli zmienię ramkę danych na Web Mercator). 3857

W przypadku EPSG (900913) dane wyjściowe są georeferencyjne, ale przesunięte o 3 komórki rastrowe na północ: 900913

Po ponownym odrzuceniu rastra za pomocą ArcGIS (eksport do WGS_1984_Web_Mercator_Auxiliary_Sphere) wynik jest prawie w porządku: Arcgis

A kiedy używam starego kodu 102113 (41001,54004) wynik jest idealny: 54004

Podsumowanie moich testów przy użyciu wszystkich kodów :

3857: far away up (missing georeference)
3785: far away up (like 3857)
3587: far away right
900913: slightly jumped up
102100: python error
102113: perfect
41001: perfect
54004: perfect
ArcGIS (web merc. aux.): good

Więc moje pytania to:

  • Dlaczego prawidłowy kod EPSG daje błędne wyniki?
  • I dlaczego stare kody działają dobrze, czy nie są przestarzałe?
  • Może moja wersja GDAL nie jest dobra lub mam błędy w kodzie Pythona?

Kod:

    yres = round(lons[1]-lons[0], 4)  # pixel size, degrees
    xres = round(lats[1]-lats[0], 4)
    ysize = len(lats)-1  # number of pixels
    xsize = len(lons)-1
    ulx = round(lons[0], 4)
    uly = round(lats[-1], 4)  # last
    driver = gdal.GetDriverByName(fileformat)
    ds = driver.Create(filename, xsize, ysize, 2, gdal.GDT_Float32)  # 2 bands
    #--- Geographic ---
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # Geographic lat/lon WGS 84
    ds.SetProjection(srs.ExportToWkt())
    gt = [ulx, xres, 0, uly, 0, -yres]  # the affine transformation coeffs (ulx, pixel, angle(skew), uly, angle, -pixel)
    ds.SetGeoTransform(gt)  # coords of top left corner of top left pixel (w-file - center of the pixel!)
    outband = ds.GetRasterBand(1)
    outband.WriteArray(data)
    outband2 = ds.GetRasterBand(2)
    outband2.WriteArray(data3)
    #--- REPROJECTION ---
    utm29 = osr.SpatialReference()
#    utm29.ImportFromEPSG(32629)  # utm 29
    utm29.ImportFromEPSG(900913)  # web mercator 3857
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(4326)
    tx = osr.CoordinateTransformation(wgs84,utm29)
    # Get the Geotransform vector
    # Work out the boundaries of the new dataset in the target projection
    (ulx29, uly29, ulz29) = tx.TransformPoint(ulx, uly)  # corner coords in utm meters
    (lrx29, lry29, lrz29) = tx.TransformPoint(ulx + xres*xsize, uly - yres*ysize )
    filenameutm = filename[0:-4] + '_web.tif'
    dest = driver.Create(filenameutm, xsize, ysize, 2, gdal.GDT_Float32)
    xres29 = round((lrx29 - ulx29)/xsize, 2) # pixel size, utm meters
    yres29 = abs(round((lry29 - uly29)/ysize, 2))
    new_gt = [ulx29, xres29, 0, uly29, 0, -yres29]
    dest.SetGeoTransform(new_gt)
    dest.SetProjection(utm29.ExportToWkt())
    gdal.ReprojectImage(ds, dest, wgs84.ExportToWkt(), utm29.ExportToWkt(), gdal.GRA_Bilinear)
    dest.GetRasterBand(1).SetNoDataValue(0.0)
    dest.GetRasterBand(2).SetNoDataValue(0.0)
    dest = None  # Flush the dataset to the disk
    ds = None  # only after the reprojected!
    print 'Image Created'
python gdal coordinate-system epsg web-mercator nadya
źródło
Może to pomóc w tym, co powiem: ponownie rzutuję raster z EPSG: 3042 na Google Mercator, pomyślałem, że to w zasadzie 3857, ale kiedy spróbuję: gdal_translate -a_srs EPSG: 3857 input.tif output.tif, wyjście jest daleko w dół (GDAL 1.11.2), na szczęście, gdy wypaczamy je za pomocą ArcGIS 10.2 i WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857 Authority: EPSG), obrazy rastrowe są we właściwym miejscu. Uważam więc, że EPSG: 3857 nie jest poprawnie obsługiwany w najnowszych wersjach GDAL.
Przedsiębiorca Web-GIS,
3
Po ponownym odrzuceniu raster nie musi już być prostokątem. Ponowne rzutowanie współrzędnych narożnych może być złym rozwiązaniem. Czy próbowałeś gdalwarp w wierszu poleceń? BTW, możesz pobrać najnowszą wersję GDAL z gisinternals.
AndreJ

Odpowiedzi:

5

Ponownie odrzuciłbym pliki za pomocą gdalwarp.

Zrobiłem to samo dla plików w EPSG: 3763, które chcę przekonwertować na EPSG: 3857. Porównałem wyniki za pomocą QGIS i Geoserver, a wygenerowane obrazy były w porządku. Ponieważ do obrazów zastosowano niewielki obrót, na granicy mogą pojawić się czarne linie (ale później linie te można uczynić przezroczystymi).

Ponieważ masz kilka tifobrazów, możesz użyć takiego skryptu , który nie zmienia żadnego istniejącego pliku, i umieszcza wygenerowane pliki w folderze o nazwie 3857:

#!/bin/bash
mkdir 3857
for file in $(ls *.tif); do
    gdalwarp -s_srs EPSG:3763 -t_srs EPSG:3857 $file 3857/$file;
    listgeo -tfw 3857/$file;
done

Jeśli chcesz również wygenerować .twfpliki, dodałem listgeo.

Ten skrypt jest przeznaczony dla systemu Linux, ale możesz napisać coś podobnego dla systemu Windows.

jgrocha
źródło