GDAL python ile başka bir ızgaraya eşleştirmek için bir ızgara nasıl projelendirilir ve yeniden örneklenir?

Question

Açıklama: Bir şekilde anahtar konu dışında kaldım: os.system veya altprocess kullanma - sadece python API'sı.

Dikey veri dönüşümleri için bir NOAA GTX ofset kılavuzunun bir bölümünü dönüştürmeye çalışıyorum ve bunu GDAL'da python ile nasıl yapacağınızı tamamen takip etmiyorum. Bir ızgara almak istiyorum (bu durumda bir Bathymetry Attributed Grid, ancak bir geotif olabilir) ve yapmak istediğim şablon olarak kullanıyorum. Bunu doğru yapabilirsem, insanların bu türden verileri kullanmasına büyük ölçüde yardımcı olacağını hissediyorum.

İşte sahip olduğum şey kesinlikle çalışmıyor. Gdalinfo'yu sonuçtaki hedef veri kümesinde (dst_ds) çalıştırdığımda, BAG kaynak kılavuzuna uymuyor.

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) 

Örnek dosyaları (ancak yapacağını, değişik çıkıntıların üst üste, ancak hiçbir iki ızgaraları):

• 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

Yapmaya çalıştığım şeye eşdeğer komut satırı:

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} 

Answer 1

Cevabınız için Jamie'ye teşekkürler.

#!/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 

Answer 2

Soruyu doğru anlıyorsam, gdalwarp ve gdal_translate alt işlemleri olarak çalıştırarak hedefinizi gerçekleştirebilirsiniz.

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) 

En zarif çözüm olmayabilir, ama bu işi olacaktır: Tam o örneğin aşağıdakileri yapın seçeneklerinizi birleştirin. Çalıştırıldıktan sonra, verilerinizi gdal kullanarak numpy'ye yükleyin ve yolunuza devam edin.