Reprojection de WGS 1984 Web Mercator (EPSG: 3857) en Python avec GDAL

17

Je reprojecte des rasters en python en utilisant GDAL. J'ai besoin de projeter plusieurs tiffs à partir des coordonnées géographiques WGS 84 vers WGS 1984 Web Mercator (Sphère auxiliaire), afin de les utiliser plus tard dans Openlayers avec OpenStreetMap et peut-être Google Maps. J'utilise Python 2.7.5 et GDAL 1.10.1 d' ici , et je transforme les coordonnées en utilisant des conseils d' ici (mon code est ci-dessous). En bref, j'ai importé osgeo.osr et utilisé ImportFromEPSG (code) et CoordinateTransformation (de, à) .

J'ai d'abord essayé EPSG (32629) qui est la zone UTM 29 et j'ai obtenu ce raster projeté (plus ou moins fin), donc le code semble être correct: utm Ensuite, j'ai utilisé EPSG (3857) parce que j'ai lu ceci et ces questions et trouvé que c'est le bon code valide récent . Mais le raster est créé sans aucune référence spatiale. Il est loin dans la trame de données WGS 84 (mais ce sera correct si je passe la trame de données à Web Mercator). 3857

Avec EPSG (900913), la sortie est géoréférencée, mais décalée d'environ 3 cellules raster vers le nord: 900913

Lorsque je reprojete le raster à l'aide d' ArcGIS (exportation dans WGS_1984_Web_Mercator_Auxiliary_Sphere), le résultat est presque correct: arcgis

Et quand j'utilise l'ancien code 102113 (41001,54004) le résultat est parfait: 54004

Le résumé de mes tests utilisant tous les 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

Mes questions sont donc:

  • Pourquoi le bon code EPSG me donne-t-il de mauvais résultats?
  • Et pourquoi les anciens codes fonctionnent bien, ne sont-ils pas obsolètes?
  • Peut-être que ma version de GDAL n'est pas bonne ou que j'ai des erreurs dans mon code python?

Le 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
la source
Cela pourrait aider ce que je vais dire, je re-projette un raster de EPSG: 3042 vers Google Mercator, je pensais être en principe le 3857, mais quand j'essaye: gdal_translate -a_srs EPSG: 3857 input.tif output.tif, la sortie est loin vers le bas (GDAL 1.11.2), heureusement lorsque vous les déformez en utilisant ArcGIS 10.2 et WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857 Authority: EPSG) les images raster sont au bon endroit. Donc, je pense qu'EPSG: 3857 n'est pas correctement géré dans les dernières versions de GDAL.
Web-GIS entrepreneur
3
Après la reprojection, le raster n'a plus besoin d'être un rectangle. Donc, reprojeter les coordonnées des coins peut être la mauvaise solution. Avez-vous essayé gdalwarp sur la ligne de commande? BTW, vous pouvez obtenir la dernière version GDAL de gisinternals.
AndreJ

Réponses:

5

Je reprojeterais les fichiers avec gdalwarp.

J'ai fait la même chose pour les fichiers EPSG: 3763 que je souhaite convertir en EPSG: 3857. J'ai comparé les résultats en utilisant QGIS et Geoserver et les images générées étaient bien. Étant donné qu'une petite rotation est appliquée aux images, vous pouvez obtenir des lignes noires sur la bordure (mais ces lignes peuvent être rendues transparentes par la suite).

Puisque vous avez plusieurs tifimages, vous pouvez utiliser un script comme celui-ci qui ne modifie aucun fichier existant et place les fichiers générés dans un dossier appelé 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

Si vous souhaitez également générer les .twffichiers, j'ai ajouté listgeo.

Ce script est pour Linux, mais vous pouvez écrire quelque chose de similaire pour Windows.

jgrocha
la source