如何投影和重采样一个网格，以匹配另一个网格与GDAL python？
如何投影和重采样一个网格，以匹配另一个网格与GDAL python？
提问于 2012-05-04 18:34:32
澄清:我忽略了关键方面:不使用os.system或子流程--只是python。
我正在尝试将NOAA GTX偏移网格的一部分转换为垂直数据转换，而不是完全遵循如何使用python在GDAL中这样做。我想使用一个网格(在本例中是一个测深属性网格，但它可能是一个geotif)，并使用它作为我想要做的模板。如果我能做到这一点，我有一种感觉，它将极大地帮助人们使用这类数据。
这是我所拥有的绝对不起作用的东西。当我在生成的目标数据集(dst_ds)上运行gdalinfo时，它与源网格包不匹配。
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)示例文件(但是任何两个网格重叠，但在不同的预测中都可以)：
- 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
这个命令行相当于我试图做的事情：
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}回答 2
Stack Overflow用户
回答已采纳
发布于 2012-05-10 16:42:04
感谢杰米的回答。
#!/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()
# 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 # FlushStack Overflow用户
发布于 2012-05-09 15:36:12
如果我正确理解这个问题，您可以通过运行gdalwarp和gdal_translate作为子进程来实现您的目标。只需组装您的选项，然后执行以下操作，例如：
import subprocess
param = ['gdalwarp',option1,option2...]
cmd = ' '.join(param)
process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout = ''.join(process.stdout.readlines())
stderr = ''.join(process.stderr.readlines())
if len(stderr) > 0:
raise IOError(stderr)这可能不是最优雅的解决方案，但它将完成工作。运行后，只需使用gdal将数据加载到numpy中，然后继续运行。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/10454316
复制相关文章
相似问题
腾讯云开发者
