Neuprojektierung von WGS 1984 Web Mercator (EPSG: 3857) in Python mit GDAL

17

Ich projiziere Raster in Python mit GDAL neu. Ich muss mehrere Tiffs von geografischen WGS 84-Koordinaten auf WGS 1984 Web Mercator (Auxiliary Sphere) projizieren, um sie später in Openlayers zusammen mit OpenStreetMap und möglicherweise Google Maps zu verwenden. Ich verwende von hier aus Python 2.7.5 und GDAL 1.10.1 und transformiere Koordinaten mithilfe von Hinweisen von hier (mein Code ist unten). Kurz gesagt, ich habe osgeo.osr importiert und ImportFromEPSG (Code) und CoordinateTransformation (from, to) verwendet .

Ich habe zuerst EPSG (32629) ausprobiert, das UTM-Zone 29 ist, und dieses projizierte Raster (mehr oder weniger gut) erhalten. Der Code scheint also korrekt zu sein: utm Dann habe ich EPSG (3857) verwendet, weil ich dieses und dieses gelesen und gefunden habe dass es der korrekte aktuelle gültige Code ist . Das Raster wird jedoch ohne räumlichen Bezug erstellt. Es ist in WGS 84-Datenrahmen weit oben (aber in Ordnung, wenn ich den Datenrahmen auf Web Mercator umstelle). 3857

Mit EPSG (900913) wird die Ausgabe georeferenziert, jedoch um 3 Rasterzellen nach Norden verschoben: 900913

Wenn ich das Raster mit ArcGIS neu projiziere (Export in WGS_1984_Web_Mercator_Auxiliary_Sphere), ist das Ergebnis nahezu in Ordnung: arcgis

Und wenn ich den alten Code 102113 (41001,54004) verwende, ist das Ergebnis perfekt: 54004

Die Zusammenfassung meiner Tests mit allen Codes :

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

Meine Fragen sind also:

  • Warum führt der richtige EPSG-Code zu falschen Ergebnissen?
  • Und warum funktionieren die alten Codes gut, sind sie nicht veraltet?
  • Vielleicht ist meine GDAL-Version nicht gut oder ich habe Fehler in meinem Python-Code?

Der Code:

    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
quelle
Es könnte helfen, was ich sagen werde, ich projiziere ein Raster von EPSG: 3042 auf das von Google Mercator, ich dachte, es wäre im Prinzip das 3857, aber wenn ich es versuche: gdal_translate -a_srs EPSG: 3857 input.tif output.tif, die Ausgabe ist weit unten (GDAL 1.11.2), glücklicherweise sind die Rasterbilder an der richtigen Stelle, wenn sie mit ArcGIS 10.2 und der WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857 Authority: EPSG) verzerrt werden. Daher glaube ich, dass EPSG: 3857 in den neuesten Versionen von GDAL nicht richtig gehandhabt wird.
Web-GIS-Unternehmer
3
Nach der Neuprojektion muss das Raster kein Rechteck mehr sein. Eine Neuprojektion der Eckkoordinaten kann daher die falsche Lösung sein. Hast du gdalwarp auf der Kommandozeile ausprobiert? Übrigens können Sie die neueste GDAL-Version von gisinternals erhalten.
AndreJ

Antworten:

5

Ich würde die Dateien mit neu projizieren gdalwarp.

Ich habe dasselbe für Dateien in EPSG: 3763 gemacht, die ich in EPSG: 3857 konvertieren möchte. Ich habe die Ergebnisse mit QGIS und Geoserver verglichen und die generierten Bilder waren in Ordnung. Da die Bilder leicht gedreht werden, kann es vorkommen, dass am Rand schwarze Linien erscheinen (diese Linien können jedoch nachträglich transparent gemacht werden).

Da Sie mehrere tifBilder haben, können Sie ein Skript wie dieses verwenden, das keine vorhandene Datei ändert und die generierten Dateien in einen Ordner mit dem Namen 3857 legt:

#!/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

Wenn Sie auch die .twfDateien generieren möchten , habe ich hinzugefügt listgeo.

Dieses Skript ist für Linux, aber Sie können etwas Ähnliches für Windows schreiben.

jgrocha
quelle