Làm thế nào để dự án và định lại mẫu lưới để khớp với lưới khác với python GDAL?

Question

Làm rõ: Tôi bằng cách nào đó đã bỏ qua khía cạnh chính: không sử dụng os.system hoặc subprocess - chỉ là python API.

Tôi đang cố chuyển đổi một phần của lưới bù đắp NOX GTX cho các phép biến đổi dữ liệu dọc và không hoàn toàn theo cách thực hiện điều này trong GDAL bằng python. Tôi muốn lấy một lưới (trong trường hợp này là Lưới phân bổ độ sâu, nhưng nó có thể là một dấu địa lý) và sử dụng nó làm mẫu mà tôi muốn làm. Nếu tôi có thể làm điều này đúng, tôi có cảm giác rằng nó sẽ giúp mọi người tận dụng loại dữ liệu này.

Đây là những gì tôi có mà chắc chắn không hoạt động. Khi tôi chạy gdalinfo trên tập dữ liệu đích (dst_ds), nó không khớp với BAG lưới nguồn.

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) 

file Ví dụ (nhưng bất kỳ hai lưới nơi họ chồng chéo lên nhau, nhưng trong dự báo khác nhau sẽ làm gì):

• 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

Dòng lệnh tương đương với những gì tôi đang cố gắng làm:

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

Nhờ Jamie cho câu trả lời.

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

Nếu tôi hiểu câu hỏi một cách chính xác, bạn có thể đạt được mục tiêu của mình bằng cách chạy gdalwarp và gdal_translate làm quy trình con. Chỉ cần tập hợp các tùy chọn của bạn, sau đó thực hiện như sau:

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) 

Đây có thể không phải là giải pháp thanh lịch nhất, nhưng nó sẽ hoàn thành công việc. Một khi nó được chạy, chỉ cần tải dữ liệu của bạn vào numpy bằng cách sử dụng gdal và thực hiện theo cách của bạn.