将shapefile与栅格对齐并将值分配给叠加层，然后作为数组返回？

Question

我的目标是将 shapefile 与栅格底图对齐，并将 1 分配给重叠的单元格，将 0 分配给不重叠的单元格，最终返回一个包含纬度、经度、时间和二进制变量 ( 1/0)。

计划如下：1) 从数组创建区域栅格，2) 栅格化多边形 shapefile，3) 将栅格化 shapefile 与基本栅格对齐，4) 重叠的像素将分配为 1，不重叠的像素将分配为 0， 5) 将栅格转换为数组。

我已经能够执行第 1 步和第 2 步（请参见下面的代码），但我在第 3 步上停留了很长时间。如何对齐两个栅格？

您可以在此处找到文件： https://www.dropbox.com/sh/pecptfepac18s2y/AADbxFkKWlLqMdiHh-ICt4UYa?dl=0

这是我用来创建 BC 平面网格作为底图的代码：

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

这是我用来栅格化多边形 shapefile 的代码

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 做所有事情，因为我无法访问或资助付费 GIS 软件。我对地理空间数据处理完全陌生……不确定我是否采取了正确的方法。任何帮助都会很棒。

Answer 1

从 rasterio 库中检出“rasterio.mask.mask”。我认为这会有所帮助。