Reproyectar 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 usarlas 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 de aquí , y transformo las coordenadas usando los consejos de < a href="http://jgomezdans.github.io/gdal_notes/reprojection.html"> aquí (mi código está debajo). En resumen, importé osgeo.osr y utilicé ImportFromEPSG (código) y CoordinateTransformation (from, to) .

Primero probé EPSG (32629) que es UTM zona 29 y obtuve este ráster proyectado (más o menos fino), por lo que el código parece ser correcto: LuegoutilicéEPSG(3857)porqueleí this y esta preguntas y encontré que es el código válido reciente correcto. Pero el raster 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).

ConEPSG(900913)lasalidaestágeorreferenciada,perosedesplazaaproximadamente3celdasrasterhaciaelnorte:

Cuando reproyecto el ráster utilizando ArcGIS (exportar en WGS_1984_Web_Mercator_Auxiliary_Sphere) el resultado es casi correcto:

Ycuandousoelcódigoanterior102113(41001,54004)elresultadoesperfecto:

El resumen de mis pruebas utilizando 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 códigos antiguos 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'
    
pregunta nadya 22.01.2014 - 00:41

2 respuestas

4

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 utilizando QGIS y Geoserver y las imágenes generadas estaban bien. Como se aplica una pequeña rotación a las imágenes, es posible que aparezcan algunas líneas negras en el borde (pero estas líneas se pueden hacer transparentes más adelante).

Ya que tiene varias imágenes tif , 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 quieres generar los archivos .twf , he agregado listgeo .

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

    
respondido por el jgrocha 26.01.2016 - 23:23
1

También iría por GDALwarp. Asegúrese de ser coherente con las "publicaciones" y las interpretaciones de "celda" de los rásteres. enlace

    
respondido por el RutgerH 02.02.2016 - 18:04

Lea otras preguntas en las etiquetas