Memproyeksikan ulang WGS 1984 Web Mercator (EPSG: 3857) dalam Python dengan GDAL

17

Saya memproyeksikan ulang raster dalam python menggunakan GDAL. Saya perlu memproyeksikan beberapa tiff dari koordinat geografis WGS 84 ke WGS 1984 Web Mercator (Auxiliary Sphere), untuk menggunakannya nanti di Openlayers bersama dengan OpenStreetMap dan mungkin Google maps. Saya menggunakan Python 2.7.5 dan GDAL 1.10.1 dari sini , dan mengubah koordinat menggunakan saran dari sini (kode saya di bawah). Singkatnya, saya mengimpor osgeo.osr dan menggunakan ImportFromEPSG (kode) dan CoordinateTransformation (dari, ke) .

Saya pertama kali mencoba EPSG (32629) yang merupakan zona UTM 29 dan mendapatkan raster yang diproyeksikan ini (kurang lebih baik), jadi kodenya tampaknya benar: utm Kemudian saya menggunakan EPSG (3857) karena saya sudah membaca ini dan ini pertanyaan dan menemukan bahwa itu adalah kode valid terbaru yang benar . Tetapi raster dibuat tanpa referensi spasial sama sekali. Itu jauh di bingkai data WGS 84 (tapi akan baik-baik saja jika saya beralih bingkai data ke Web Mercator). 3857

Dengan EPSG (900913) output di-georeferensi, tetapi bergeser sekitar 3 sel raster ke utara: 900913

Ketika saya memproyeksikan ulang raster menggunakan ArcGIS (ekspor di WGS_1984_Web_Mercator_Auxiliary_Sphere) hasilnya hampir baik: arcgis

Dan ketika saya menggunakan kode lama 102113 (41001.54004) hasilnya sempurna: 54004

Ringkasan tes saya menggunakan semua kode :

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

Jadi pertanyaan saya adalah:

  • Mengapa kode EPSG yang benar memberi saya hasil yang salah?
  • Dan mengapa kode lama berfungsi dengan baik, bukankah itu sudah usang?
  • Mungkin versi GDAL saya tidak baik atau saya memiliki kesalahan dalam kode python saya?

Kode:

    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
sumber
Mungkin membantu apa yang akan saya katakan, saya memproyeksikan ulang raster dari EPSG: 3042 ke Google Mercator, saya pikir pada prinsipnya 3857, tetapi ketika saya mencoba: gdal_translate -a_srs EPSG: 3857 input.tif output.tif, outputnya jauh di bawah (GDAL 1.11.2), untungnya ketika melengkungkan mereka menggunakan ArcGIS 10.2 dan WGS_1984_Web_Mercator_Auxiliary_Sphere (WKID: 3857 Authority: EPSG) gambar raster berada di tempat yang tepat. Jadi, saya percaya EPSG: 3857 tidak ditangani dengan benar dalam versi terbaru GDAL.
Pengusaha Web-GIS
3
Setelah proyeksi ulang, raster tidak harus persegi panjang lagi. Jadi memproyeksikan kembali koordinat sudut mungkin merupakan solusi yang salah. Sudahkah Anda mencoba gdalwarp di baris perintah? BTW Anda bisa mendapatkan versi GDAL terbaru dari gisinternals.
AndreJ

Jawaban:

5

Saya akan memproyeksikan ulang file dengan gdalwarp.

Saya telah melakukan hal yang sama untuk file dalam EPSG: 3763 yang ingin saya konversi ke EPSG: 3857. Saya membandingkan hasil menggunakan QGIS dan Geoserver dan gambar yang dihasilkan baik-baik saja. Karena rotasi kecil diterapkan pada gambar, Anda mungkin mendapatkan beberapa garis hitam di perbatasan (tetapi garis-garis ini dapat dibuat transparan setelahnya).

Karena Anda memiliki beberapa tifgambar, Anda dapat menggunakan skrip seperti ini yang tidak mengubah file apa pun yang ada, dan menempatkan file yang dihasilkan dalam folder bernama 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

Jika Anda juga ingin membuat .twffile, saya telah menambahkan listgeo.

Skrip ini untuk Linux, tetapi Anda dapat menulis sesuatu yang serupa untuk Windows.

jgrocha
sumber