Как спроецировать и передискретизировать сетку, чтобы она соответствовала другой сетке с помощью GDAL python?

Пояснение :Я почему-то упустил ключевой аспект :не используя os.system или подпроцесс -только python API.

Я пытаюсь преобразовать часть сетки смещения NOAA GTX для преобразования вертикальных датумов и не совсем понимаю, как это сделать в GDAL с помощью python. Я хотел бы взять сетку (, в данном случае сетку с атрибутами батиметрии, но это может быть геотиф )и использовать его в качестве шаблона, который я хотел бы сделать. Если я смогу сделать это правильно, у меня есть ощущение, что это очень поможет людям использовать этот тип данных.

Вот что у меня определенно не работает. Когда я запускаю gdalinfo для результирующего целевого набора данных (dst _ds ), он не соответствует исходной сетке BAG.

from osgeo import gdal, osr

bag = gdal.Open(bag_filename)
gtx = gdal.Open(gtx_filename)

bag_srs = osr.SpatialReference()
bag_srs.ImportFromWkt(bag.GetProjection())

vrt = gdal.AutoCreateWarpedVRT(gtx, None, bag_srs.ExportToWkt(), gdal.GRA_Bilinear,  0.125)

dst_ds = gdal.GetDriverByName('GTiff').Create(out_filename, bag.RasterXSize, bag.RasterYSize,
                                            1, gdalconst.GDT_Float32)
dst_ds.SetProjection(bag_srs.ExportToWkt())
dst_ds.SetGeoTransform(vrt.GetGeoTransform())

def warp_progress(pct, message, user_data):
  return 1

gdal.ReprojectImage(gtx, dst_ds, None, None, gdal.GRA_NearestNeighbour, 0, 0.125, warp_progress, None)

Примеры файлов (, но подойдут любые две сетки, в которых они перекрываются, но находятся в разных проекциях):

Командная строка, эквивалентная тому, что я пытаюсь сделать:

gdalwarp -tr 2 -2 -te 369179 4773093 372861 4775259 -of VRT -t_srs EPSG:2960 \
     MENHMAgome01_8301/mllw.gtx  mllw-2960-crop-resample.vrt
gdal_translate mllw-2960-crop-resample.{vrt,tif}

12
задан Wraf 2 October 2013 в 12:38
поделиться