Reprojetando o WGS 1984 Web Mercator (EPSG: 3857) em Python com GDAL

17

Estou reprojetando rasters em python usando GDAL. Preciso projetar vários tiffs, desde as coordenadas geográficas WGS 84 até o WGS 1984 Web Mercator (Auxiliary Sphere), para usá-los posteriormente em Openlayers, juntamente com o OpenStreetMap e talvez com o Google maps. Estou usando Python 2.7.5 e GDAL 1.10.1 a partir daqui e transformando coordenadas usando conselhos daqui (meu código está abaixo). Em resumo, importei o osgeo.osr e usei ImportFromEPSG (código) e CoordateTransformation (de, para) .

Tentei primeiro o EPSG (32629), que é a zona 29 da UTM, e obtive essa varredura projetada (mais ou menos fina); portanto, o código parece estar correto: utm Então usei o EPSG (3857) porque li estas e estas perguntas e encontrei que é o código válido recente correto . Mas a varredura é criada sem nenhuma referência espacial. Está muito distante no quadro de dados WGS 84 (mas ficará bem se eu trocar o quadro de dados para Web Mercator). 3857

Com o EPSG (900913), a saída é georreferenciada, mas deslocou cerca de 3 células raster para o norte: 900913

Quando reprojeto a varredura usando o ArcGIS (exportação no WGS_1984_Web_Mercator_Auxiliary_Sphere), o resultado é quase bom: arcgis

E quando eu uso o código antigo 102113 (41001,54004), o resultado é perfeito: 54004

O resumo dos meus testes usando todos os 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

Então, minhas perguntas são:

  • Por que o código EPSG correto me fornece resultados errados?
  • E por que os códigos antigos funcionam bem, eles não são preteridos?
  • Talvez minha versão GDAL não seja boa ou eu tenha erros no meu código python?

O 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
fonte
Pode ajudar o que eu vou dizer, estou reprojetando uma varredura do EPSG: 3042 para a do Google Mercator, que eu pensava ser o 3857, mas quando tento: gdal_translate -a_srs EPSG: 3857 input.tif output.tif, a saída está muito distante (GDAL 1.11.2), felizmente quando os deformamos usando o ArcGIS 10.2 e o WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857 Authority: EPSG) as imagens rasterizadas estão no lugar certo. Portanto, acredito que o EPSG: 3857 não é tratado adequadamente nas versões mais recentes do GDAL.
Empresário Web-GIS
3
Após a reprojeção, a varredura não precisa mais ser um retângulo. Portanto, reprojetar as coordenadas dos cantos pode ser a solução errada. Você já tentou o gdalwarp na linha de comando? BTW, você pode obter a versão mais recente do GDAL na gisinternals.
11556 AndreJ

Respostas:

5

Eu reprojetaria os arquivos com gdalwarp.

Fiz o mesmo para arquivos no EPSG: 3763 que desejo converter para EPSG: 3857. Eu comparei os resultados usando QGIS e Geoserver e as imagens geradas foram boas. Como uma pequena rotação é aplicada às imagens, você pode obter algumas linhas pretas na borda (mas essas linhas podem ficar transparentes depois).

Como você possui várias tifimagens, é possível usar um script como este que não altera nenhum arquivo existente e coloca os arquivos gerados em uma pasta chamada 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

Se você também deseja gerar os .twfarquivos, eu adicionei listgeo.

Este script é para Linux, mas você pode escrever algo semelhante para o Windows.

jgrocha
fonte