Обрезка узлов из пространственного изображения - PullRequest
Новая
       18

Обрезка узлов из пространственного изображения

1 голос
/ 08 марта 2019

У меня есть несколько растров, где информация о пикселях - это данные высот. Однако в изображении много значений ноданных (нулей), выходящих за пределы объема данных, которые мне нужны, поэтому размеры файлов неоправданно велики. Я попытался обрезать изображение, используя сочетание gdal и numpy, но я заметил небольшое изменение размера субпикселя в изображении. Это ошибка округления, и можно ли что-нибудь сделать, чтобы решить эту проблему?

В некоторых местах он не очень питоничен, но, кажется, работает нормально, за исключением небольшой смены

Я также открыт для любых альтернативных способов сделать это с помощью python / numpy / gdal ...

Мой код: 

from osgeo import gdal
import numpy as np
import os
import subprocess

np.set_printoptions(threshold=np.nan)

inDir = '/path/to/files'
outDir = '/path/for/outfiles'

files = [file for file in os.listdir(inDir) if file.endswith('.tif')]

for file in files:
    outfile= os.path.join(outDir, file)
    if os.path.isfile(outfile):
        print('File exists')
        pass
    else:
        try:
            print(file)

            # open raster and read in the data
            dataset = gdal.Open(file)
            band = dataset.GetRasterBand(1)
            Cols = dataset.RasterXSize
            Rows = dataset.RasterYSize
            print('Reading Radar data...')
            data = band.ReadAsArray(0, 0, Cols, Rows).astype(np.float)
            # Read in the projection onfo
            ulx, xres, xskew, uly, yskew, yres  = dataset.GetGeoTransform()

            # define the corner co-ords
            ulx = ulx
            uly = uly
            lrx = ulx + (dataset.RasterXSize * xres)
            lry = uly + (dataset.RasterYSize * yres)


            # Make two arrays with the co-ords for each pixel
            XP = np.arange(ulx, lrx, xres, dtype=np.float64)
            YP = np.arange(uly, lry, yres, dtype=np.float64)            

            # Make an array the same shape as data populated with the X dimensions
            print('Building X dimensions...')
            Xarr = np.zeros_like(data)
            Xarr[...] = XP

            # Make an array the same shape as data populated with the Y dimensions
            print('Building Y dimensions...')
            Yarr = np.zeros_like(data)
            Yarr = Yarr.transpose()
            Yarr[...] = YP
            Yarr = Yarr.transpose()

            # Mask out the X and Y co-ords arrays based on presence of valid data
            print('Subsetting data...')
            SubXarr = Xarr[data > 0]
            SubYarr = Yarr[data > 0]
            SubData = data[data > 0]

            #Get the bounding co-ords based on presence of valid data
            # note yres is -ve
            print('Finding corner coords...')
            Xmin = np.amin(SubXarr)
            Xmax = np.amax(SubXarr) + abs(xres)
            Ymin = np.amin(SubYarr) - abs(yres)
            Ymax = np.amax(SubYarr)

            print(Xmin, Ymin)
            print(Xmax, Ymax)

            # Subset the input image based on the bounding co-ords
            outfile = os.path.join(outDir, file)
            print('Cropping Data to ' + str(outfile.split('/')[-1]))

            subprocess.call('gdalwarp -of GTiff -r cubic -tr ' + str(xres) + ' ' + str(yres) + ' -te ' + str(Xmin) + ' ' + str(Ymin) + ' ' + str(Xmax) + ' ' + str(Ymax) + ' -tr ' + str(xres) + ' ' + str(yres) + ' ' + str(file) + ' ' + str(outfile),shell=True)

        except Exception as e:
            print(e)
            print('Could not process: ', str(outfile.split('/')[-1]))
...