Reproyección de WGS 1984 Web Mercator (EPSG: 3857) en Python con GDAL

17

Estoy reproyectando rásteres en python usando GDAL. Necesito proyectar varios tiffs desde las coordenadas geográficas WGS 84 a WGS 1984 Web Mercator (esfera auxiliar), para poder usarlos más tarde en Openlayers junto con OpenStreetMap y tal vez los mapas de Google. Estoy usando Python 2.7.5 y GDAL 1.10.1 desde aquí , y transformando coordenadas usando consejos desde aquí (mi código está abajo). En resumen, importé osgeo.osr y usé ImportFromEPSG (código) y CoordinateTransformation (de, a) .

Primero probé EPSG (32629), que es la zona 29 de UTM y obtuve este ráster proyectado (más o menos fino), por lo que el código parece ser correcto: utm luego usé EPSG (3857) porque leí esto y estas preguntas y encontré que es el código válido reciente correcto . Pero el ráster se crea sin ninguna referencia espacial. Está muy lejos en el marco de datos WGS 84 (pero estará bien si cambio el marco de datos a Web Mercator). 3857

Con EPSG (900913), la salida está georreferenciada, pero se desplaza unas 3 celdas ráster hacia el norte: 900913

Cuando reproyecto el ráster usando ArcGIS (exportar en WGS_1984_Web_Mercator_Auxiliary_Sphere) el resultado es casi bueno: arcgis

Y cuando uso el antiguo código 102113 (41001,54004) el resultado es perfecto: 54004

El resumen de mis pruebas con todos los códigos :

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

Entonces mis preguntas son:

  • ¿Por qué el código EPSG correcto me da resultados incorrectos?
  • ¿Y por qué los viejos códigos funcionan bien, no están en desuso?
  • ¿Quizás mi versión GDAL no es buena o tengo errores en mi código de Python?

El código:

    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
fuente
Podría ayudar lo que voy a decir, estoy re-proyectando un ráster de EPSG: 3042 al de Google Mercator, pensé que en principio era el 3857, pero cuando intento: gdal_translate -a_srs EPSG: 3857 input.tif output.tif, la salida está muy lejos (GDAL 1.11.2), afortunadamente cuando los deforma con ArcGIS 10.2 y WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857 Authority: EPSG) las imágenes ráster están en el lugar correcto. Por lo tanto, creo que EPSG: 3857 no se maneja adecuadamente en las últimas versiones de GDAL.
Web-GIS emprendedor
3
Después de la reproyección, el ráster ya no tiene que ser un rectángulo. Entonces, reproyectar las coordenadas de las esquinas puede ser la solución incorrecta. ¿Has probado gdalwarp en la línea de comando? Por cierto, puede obtener la última versión GDAL de gisinternals.
AndreJ

Respuestas:

5

Reproyectaría los archivos con gdalwarp.

He hecho lo mismo para los archivos en EPSG: 3763 que quiero convertir a EPSG: 3857. Comparé los resultados usando QGIS y Geoserver y las imágenes generadas estaban bien. Dado que se aplica una pequeña rotación a las imágenes, puede obtener algunas líneas negras en el borde (pero estas líneas pueden hacerse transparentes después).

Como tiene varias tifimágenes, puede usar un script como este que no cambia ningún archivo existente y coloca los archivos generados en una carpeta llamada 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 también desea generar los .twfarchivos, he agregado listgeo.

Este script es para Linux, pero puede escribir algo similar para Windows.

jgrocha
fuente