Я репроектирую растры в Python, используя GDAL. Мне нужно спроецировать несколько tiff из географических координат WGS 84 в WGS 1984 Web Mercator (Вспомогательная сфера), чтобы потом использовать их в Openlayers вместе с OpenStreetMap и, возможно, картами Google. Я использую Python 2.7.5 и GDAL 1.10.1 отсюда , и преобразовываю координаты, используя советы отсюда (мой код ниже). Короче говоря, я импортировал osgeo.osr и использовал ImportFromEPSG (код) и CoordinateTransformation (от, до) .
Сначала я попробовал EPSG (32629), который является зоной UTM 29, и получил этот проецируемый растр (более или менее хорошо), поэтому код кажется правильным: Потом я использовал EPSG (3857), потому что я прочитал этот и этот вопросы и нашел что это правильный недавний действительный код . Но растр создается без какой-либо пространственной привязки. Он находится далеко во фрейме данных WGS 84 (но все будет в порядке, если я переключу фрейм данных на Web Mercator).
В EPSG (900913) выходные данные имеют географическую привязку, но смещены примерно на 3 растровых ячейки на север:
Когда я перепроектирую растр с помощью ArcGIS (экспорт в WGS_1984_Web_Mercator_Auxili_Sphere), результат почти идеален:
И когда я использую старый код 102113 (41001,54004), результат идеален:
Резюме моих тестов с использованием всех кодов :
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
Итак, мои вопросы:
- Почему правильный код EPSG дает мне неправильные результаты?
- И почему старые коды работают нормально, разве они не устарели?
- Может быть, моя версия GDAL не подходит или у меня есть ошибки в моем коде Python?
Код:
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'
Ответы:
Я бы перепроектировал файлы с
gdalwarp
.Я сделал то же самое для файлов в EPSG: 3763, которые я хочу преобразовать в EPSG: 3857. Я сравнил результаты, используя QGIS и Geoserver, и сгенерированные изображения были в порядке. Поскольку к изображениям применяется небольшое вращение, на границе могут появиться черные линии (но впоследствии эти линии можно сделать прозрачными).
Поскольку у вас есть несколько
tif
изображений, вы можете использовать такой скрипт , который не изменяет ни один из существующих файлов и помещает сгенерированные файлы в папку с именем 3857:Если вы также хотите сгенерировать
.twf
файлы, я добавилlistgeo
.Этот скрипт для Linux, но вы можете написать что-то похожее для Windows.
источник
Пойдет и на GDALwarp. Убедитесь, что вы согласны с интерпретациями «проводок» и «ячеек» растров. http://www.remotesensing.org/geotiff/spec/geotiff2.5.html
источник