Ошибка при извлечении растра по маске с помощью osgeo.gdal.Translate - PullRequest
Новая
       49

Ошибка при извлечении растра по маске с помощью osgeo.gdal.Translate

0 голосов
/ 17 апреля 2020

Я использую gdal.Translate, чтобы извлечь растр по маске. Обе растровые данные имеют одинаковый CRS и одинаковый размер пикселя. Вот фрагмент кода. 

def by_mask(in_dir, in_fname, mask_dir, mask_fname, out_dir, out_fname):    
    mask_ds = gdal.Open(os.path.join(mask_dir, mask_fname), GA_ReadOnly)
    mask_proj = mask_ds.GetProjection()
    mask_geotrans = mask_ds.GetGeoTransform()
    ulx = mask_geotrans[0]
    uly = mask_geotrans[3]
    lrx = ulx + mask_geotrans[1] * mask_ds.RasterXSize
    lry = uly + mask_geotrans[5] * mask_ds.RasterYSize

    ds = gdal.Open(os.path.join(in_dir, in_fname), GA_ReadOnly)
    gdal.Translate(destName=os.path.join(out_dir, out_fname),
                   srcDS=ds,
                   options=gdal.TranslateOptions(
                       format='GTiff',
                       projWin=[ulx, uly, lrx, lry],
                       options=['-co', 'COMPRESS=DEFLATE'],
                       xRes=res_info(in_dir=mask_dir, in_fname=mask_fname)[0],
                       yRes=res_info(in_dir=mask_dir, in_fname=mask_fname)[1]
                   ))

И я получил сообщение об ошибке ERROR 1: Error: Computed -srcwin 11577 6424 1236 -873 has negative width and/or height. Есть предложения?

1 Ответ

0 голосов
/ 26 апреля 2020

Вот мое решение. Любые предложения или вопросы приветствуются. 

from osgeo import gdal
from gdalconst import GA_ReadOnly
import os
import osr


def by_mask(in_dir, in_fname, mask_dir, mask_fname, out_dir, out_fname,
            mask_rev_lst):
    """
    extract raster data with another mask raster data, working well with
    inconsistent extents
    :param in_dir: directory of target raster data
    :type in_dir: str
    :param in_fname: filename of target raster data
    :type in_fname: str
    :param mask_dir: directory of mask raster data
    :type mask_dir: str
    :param mask_fname: filename of mask raster data
    :type mask_fname: str
    :param out_dir: directory of raster data to be saved
    :type out_dir: str
    :param out_fname: filename of raster data to be saved
    :type out_fname: str
    :param mask_rev_lst: values to masked
    :type mask_rev_lst: list
    :return:
    :rtype:
    """
    src_f = os.path.join(in_dir, in_fname)
    mask_f = os.path.join(mask_dir, mask_fname)
    dest_f = os.path.join(out_dir, out_fname)

    src_ds = gdal.Open(src_f)
    src_nodata = src_ds.GetRasterBand(1).GetNoDataValue()
    exists_nodata = False if src_nodata is None else True

    mask_ds = gdal.Open(mask_f)
    mask_arr = mask_ds.ReadAsArray()
    mask_geotrans = mask_ds.GetGeoTransform()
    mask_proj = mask_ds.GetProjection()

    mem_driver = gdal.GetDriverByName('Mem')
    mem_ds = mem_driver.Create('', mask_ds.RasterXSize, mask_ds.RasterYSize,
                               1, gdal.GDT_Float32)
    mem_ds.SetProjection(mask_proj)
    mem_ds.SetGeoTransform(mask_geotrans)
    mem_srs = osr.SpatialReference()
    mem_srs.ImportFromWkt(mask_proj)
    if src_nodata is not None:
        mem_ds.GetRasterBand(1).SetNoDataValue(src_nodata)

    gdal.Warp(destNameOrDestDS=mem_ds,
              srcDSOrSrcDSTab=src_ds,
              options=gdal.WarpOptions(
                  resampleAlg=gdal.gdalconst.GRA_NearestNeighbour,
                  srcNodata=src_nodata,
                  dstNodata=src_nodata
              ))
    mem_arr = mem_ds.ReadAsArray()
    del mem_ds
    mem_arr_lst = []
    for el in mask_rev_lst:
        temp_arr = mem_arr.copy()
        temp_arr[~(mask_arr == el)] = 1e29
        if exists_nodata:
            mem_arr_lst.append(temp_arr)
        else:
            mem_arr_lst.append(temp_arr[~(temp_arr == 1e29)])


    dest_driver = gdal.GetDriverByName('GTiff')
    dest_ds = dest_driver.Create(dest_f,
                                 mask_ds.RasterXSize, mask_ds.RasterYSize,
                                 1, gdal.GDT_Float32,
                                 options=['COMPRESS=DEFLATE'])
    dest_ds.SetGeoTransform(mask_geotrans)
    dest_ds.SetProjection(mask_proj)
    dest_ds.GetRasterBand(1).WriteArray(*mem_arr_lst)
    if src_nodata is not None:
        dest_ds.GetRasterBand(1).SetNoDataValue(1e29)

    dest_ds.FlushCache()
...