将shapefile对齐到光栅并将值分配给覆盖，然后作为数组返回？

Question

我的目标是将shapefile与光栅底图对齐，将1指定给重叠的单元格，将0指定给不重叠的单元格，最终返回一个包含lat、lon、time和二进制变量（1/0）的数组。

计划如下：1）从阵列创建区域光栅，2）光栅化多边形形状文件，3）将光栅化形状文件与基础光栅对齐，4）重叠的像素将被指定为1，不重叠的像素将被指定为0，5）将光栅转换为阵列。

我已经完成了第一步

你可以在这里找到这些文件：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

下面是我用来栅格化多边形形状文件的代码

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图书馆的“面具”。我想这会有帮助的。