如何使用GDAL python对网格进行投影和重新采样以使其与另一个网格匹配? [英] How to project and resample a grid to match another grid with GDAL python?

问题描述

说明:我不知何故忽略了关键方面:不使用os.system或子进程-仅使用python API.

Clarification: I somehow left out the key aspect: not using os.system or subprocess - just the python API.

我正在尝试转换NOAA GTX偏移网格的一部分以进行垂直基准面转换，而不是完全遵循在GDAL中使用python进行此操作的方法.我想要一个网格(在本例中为测深归因网格"，但它可能是一个geotif)，并将其用作我想要做的模板.如果我能做到这一点，我有一种感觉，它将极大地帮助人们使用这种类型的数据.

I'm trying to convert a section of a NOAA GTX offset grid for vertical datum transformations and not totally following how to do this in GDAL with python. I'd like to take a grid (in this case a Bathymetry Attributed Grid, but it could be a geotif) and use it as the template that I'd like to do to. If I can do this right, I have a feeling that it will greatly help people make use of this type of data.

这是我所拥有的绝对无法正常工作的东西.在生成的目标数据集(dst_ds)上运行gdalinfo时，它与源网格BAG不匹配.

Here is what I have that is definitely not working. When I run gdalinfo on the resulting destination dataset (dst_ds), it does not match the source grid 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)

示例文件(但重叠的两个网格但投影位置不同的两个网格都可以):

Example files (but any two grids where they overlap, but are in different projections would do):

• http://surveys.ngdc.noaa .gov/mgg/NOS/coast/F00001-F02000/F00574/BAG/ F00574_MB_2m_MLLW_2of3.bag
• http://vdatum.noaa.gov/download/data/VDatum_National.zip MENHMAgome01_8301/mllw.gtx

• http://surveys.ngdc.noaa.gov/mgg/NOS/coast/F00001-F02000/F00574/BAG/ F00574_MB_2m_MLLW_2of3.bag
• http://vdatum.noaa.gov/download/data/VDatum_National.zip MENHMAgome01_8301/mllw.gtx

等效于我要执行的命令行:

The command line equivalent to what I'm trying to do:

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}

推荐答案

感谢杰米的回答.

#!/usr/bin/env python

from osgeo import gdal, gdalconst

# Source
src_filename = 'MENHMAgome01_8301/mllw.gtx'
src = gdal.Open(src_filename, gdalconst.GA_ReadOnly)
src_proj = src.GetProjection()
src_geotrans = src.GetGeoTransform()

# We want a section of source that matches this:
match_filename = 'F00574_MB_2m_MLLW_2of3.bag'
match_ds = gdal.Open(match_filename, gdalconst.GA_ReadOnly)
match_proj = match_ds.GetProjection()
match_geotrans = match_ds.GetGeoTransform()
wide = match_ds.RasterXSize
high = match_ds.RasterYSize

# Output / destination
dst_filename = 'F00574_MB_2m_MLLW_2of3_mllw_offset.tif'
dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Float32)
dst.SetGeoTransform( match_geotrans )
dst.SetProjection( match_proj)

# Do the work
gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)

del dst # Flush

这篇关于如何使用GDAL python对网格进行投影和重新采样以使其与另一个网格匹配?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！