Выровнять шейп-файл по растру и присвоить значения наложению, а затем вернуть в виде массива? - PullRequest
Новая
       133

Выровнять шейп-файл по растру и присвоить значения наложению, а затем вернуть в виде массива?

0 голосов
/ 14 июля 2020

Моя цель - выровнять шейп-файл с растровой базовой картой и присвоить 1 перекрывающимся ячейкам и 0 тем, которые не перекрываются, в конечном итоге возвращая массив, содержащий широту, долготу, время и двоичную переменную ( 1/0).

Вот план: 1) создать растр области из массива, 2) растрировать шейп-файлы многоугольника, 3) выровнять растеризованные шейп-файлы с базовым растром, 4) пикселям, которые перекрывают это не будет 0, 5) преобразовать растры в массив.

Я смог выполнить шаги 1 и 2 (см. код ниже), но я застрял на шаге 3 надолго время. Как выровнять два растра?

Вы можете найти файлы здесь: https://www.dropbox.com/sh/pecptfepac18s2y/AADbxFkKWlLqMdiHh-ICt4UYa?dl=0

Вот код, который я использовал для создания плоской сетки B C как базовая карта: 

import gdal, osr
import numpy as np

#define parameters
#units = km
grid_size = 5
BC_width  = 700
BC_length = 1800

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    
    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()



def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    reversed_arr = array[::-1] # reverse array so the tif looks like the array
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster



if __name__ == "__main__":
    array = np.zeros([int(BC_length/grid_size),int(BC_width/grid_size)])  #140x360
    for i in range(1,100):
        array[i] = 100
    rasterOrigin = (-139.72938, 47.655534) #lower left corner of raster
    newRasterfn = '/temp/test.tif'
    cols = array.shape[1] #shape of an array (aka # of elements in each dimension)
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    pixelWidth = 5
    pixelHeight = 5

Вот код, который я использовал для растеризации шейп-файлов полигонов 

import ogr, gdal, osr

output_raster = '/testdata/poly.tif'

shapefile = "/testdata/20180808.shp"

def main(shapefile):

    #making the shapefile as an object.
    input_shp = ogr.Open(shapefile)
    
    #getting layer information of shapefile.
    shp_layer = input_shp.GetLayer()
    
    #pixel_size determines the size of the new raster.
    #pixel_size is proportional to size of shapefile.
    pixel_size = 0.1
    
    #get extent values to set size of output raster.
    x_min, x_max, y_min, y_max = shp_layer.GetExtent()
    
    #calculate size/resolution of the raster.
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)
    
    #get GeoTiff driver by 
    image_type = 'GTiff'
    driver = gdal.GetDriverByName(image_type)
    #passing the filename, x and y direction resolution, no. of bands, new raster.
    new_raster = driver.Create(output_raster, x_res, y_res, 1, gdal.GDT_Byte)
    
    #transforms between pixel raster space to projection coordinate space.
    new_raster.SetGeoTransform((x_min, pixel_size, 0, y_min, 0, pixel_size))
    
    #get required raster band.
    band = new_raster.GetRasterBand(1)
    
    #assign no data value to empty cells.
    no_data_value = -9999
    band.SetNoDataValue(no_data_value)
    band.FlushCache()
    
    #main conversion method
    gdal.RasterizeLayer(new_raster, [1], shp_layer, burn_values=[255])
    
    #adding a spatial reference
    new_rasterSRS = osr.SpatialReference()
    new_rasterSRS.ImportFromEPSG(4326)
    new_raster.SetProjection(new_rasterSRS.ExportToWkt())
    return output_raster

Я делаю все в Python, так как у меня нет доступа или финансирования платное ПО ГИС. Я новичок в обработке геопространственных данных ... не уверен, что использую правильный подход. Любая помощь была бы потрясающей.

...