Python osgeo.gdal.RasterizeLayer() Examples

The following are code examples for showing how to use osgeo.gdal.RasterizeLayer(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 6 votes vote down vote up
def lyr2mask(lyr, r_ds):
    #Create memory raster dataset and fill with 0s
    m_ds = gdal.GetDriverByName('MEM').CreateCopy('', r_ds, 1)
    b = m_ds.GetRasterBand(1)
    b.Fill(0)
    b.SetNoDataValue(0)
    #Not sure if gdal.RasterizeLayer can handle srs difference
    #r_srs = get_ds_srs(m_ds)
    #lyr_srs = lyr.GetSpatialReference()
    #if not r_srs.IsSame(lyr_srs):
    #   lyr = lyr_proj(lyr, r_srs)
    #Rasterize with values of 1
    gdal.RasterizeLayer(m_ds, [1], lyr, burn_values=[1])
    a = b.ReadAsArray()
    mask = a.astype('Bool')
    m_ds = None
    return ~mask

#Rasterize ogr geometry to binary mask (creates dummy layer) 
Example 2
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 5 votes vote down vote up
def rasterize(self, griddef=None, layer=0, invert=False, asVar=False, ldebug=False):
    ''' "burn" shapefile on a 2D raster; returns a 2D boolean array '''
    if griddef.__class__.__name__ != GridDefinition.__name__: raise TypeError 
    #if not isinstance(griddef,GridDefinition): raise TypeError # this is always False. probably due to pickling
    if not isinstance(invert,(bool,np.bool)): raise TypeError
    # fill values
    if invert: inside, outside = 0,1
    else: inside, outside = 1,0
    shp_lyr = self.getLayer(layer) # get shape layer
    # create raster to burn shape onto
    if ldebug: print(' - creating raster')
    msk_ds = ramdrv.Create(self.name, griddef.size[0], griddef.size[1], 1, gdal.GDT_Byte)
    # N.B.: this is a special case: only one band (1) and always boolean (gdal.GDT_Byte)
    # set projection parameters
    msk_ds.SetGeoTransform(griddef.geotransform)  # does the order matter?
    msk_ds.SetProjection(griddef.projection.ExportToWkt())  # is .ExportToWkt() necessary?
    # initialize raster band        
    msk_rst = msk_ds.GetRasterBand(1) # only one anyway...
    msk_rst.Fill(outside); msk_rst.SetNoDataValue(outside) # fill with zeros
    # burn shape layer onto raster band
    if ldebug: print(' - burning layer to raster')
    err = gdal.RasterizeLayer(msk_ds, [1], shp_lyr, burn_values = [inside]) # None, None, [1] # burn_value = 1
    # use argument ['ALL_TOUCHED=TRUE'] like so: None, None, [1], ['ALL_TOUCHED=TRUE']
    if err != 0: raise GDALError('ERROR CODE %i'%err)
    #msk_ds.FlushCash()  
    # retrieve mask array from raster band
    if ldebug: print(' - retrieving mask')
    mask = msk_ds.GetRasterBand(1).ReadAsArray()
    # convert to Variable object, is desired
    if asVar: 
      mask = Variable(name=self.name, units='mask', axes=(griddef.ylat,griddef.xlon), data=mask, 
                      dtype=np.bool, mask=None, fillValue=outside, atts=None, plot=None)
      mask = addGDALtoVar(mask, griddef=griddef,) # add GDAL to mask       
    # return mask array
    return mask  

# a container class that also acts as a shape for the outline 
Example 3
Project: coded   Author: bullocke   File: classify.py    MIT License 5 votes vote down vote up
def create_mask_from_vector(vector_data_path, cols, rows, geo_transform, 
			    projection, target_value=1,
                            output_fname='', dataset_format='MEM'):

    """
    Rasterize the given vector (wrapper for gdal.RasterizeLayer). 
    Return a gdal.Dataset.
    :param vector_data_path: Path to a shapefile
    :param cols: Number of columns of the result
    :param rows: Number of rows of the result
    :param geo_transform: Returned value of gdal.Dataset.GetGeoTransform 
	(coefficients for transforming between pixel/line (P,L) raster space,
	 and projection coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by 
	gdal.Dataset.GetProjectionRef)
    :param target_value: Pixel value for the pixels. Must be a valid 
	gdal.GDT_UInt16 value.
    :param output_fname: If the dataset_format is GeoTIFF, this is the output 
	file name
    :param dataset_format: The gdal.Dataset driver name. [default: MEM]
    """

    driver = ogr.GetDriverByName('ESRI Shapefile')
    data_source = driver.Open(vector_data_path, 0)
    if data_source is None:
        report_and_exit("File read failed: %s", vector_data_path)
    layer = data_source.GetLayer(0)
    driver = gdal.GetDriverByName(dataset_format)
    target_ds = driver.Create(output_fname, cols, rows, 1, gdal.GDT_UInt16)
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(projection)
    gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[target_value])
    return target_ds 
Example 4
Project: ALCD   Author: CNES   File: masks_preprocessing.py    GNU General Public License v3.0 5 votes vote down vote up
def rasterize_shp(input_shp, out_tif, reference_tif):
    '''
    From a shapefile, rasterize it
    '''
    gdalformat = 'GTiff'
    datatype = gdal.GDT_Byte
    burnVal = 0  # value for the output image pixels

    # Get projection info from reference image
    image = gdal.Open(reference_tif, gdal.GA_ReadOnly)
    print(image)

    # Open Shapefile
    shapefile = ogr.Open(input_shp)
    shapefile_layer = shapefile.GetLayer()

    # Rasterise
    output = gdal.GetDriverByName(gdalformat).Create(
        out_tif, image.RasterXSize, image.RasterYSize, 1, datatype, options=['COMPRESS=DEFLATE'])
    output.SetProjection(image.GetProjectionRef())
    output.SetGeoTransform(image.GetGeoTransform())

    # Write data to band 1
    band = output.GetRasterBand(1)
    band.SetNoDataValue(1)
    gdal.RasterizeLayer(output, [1], shapefile_layer, burn_values=[burnVal])

    # Close datasets
    band = None
    output = None
    image = None
    shapefile = None 
Example 5
Project: PyGEM   Author: drounce   File: shean_mb_parallel.py    MIT License 5 votes vote down vote up
def get_date_a(ds, date_shp_lyr, glac_geom_mask, datefield):
    date_r_ds = iolib.mem_drv.CreateCopy('', ds)
    #Shapefile order should be sorted by time, but might want to think about sorting here
    #Can automatically search for datefield
    gdal.RasterizeLayer(date_r_ds, [1], date_shp_lyr, options=["ATTRIBUTE=%s" % datefield])
    date_a = np.ma.array(iolib.ds_getma(date_r_ds), mask=glac_geom_mask)
    #Note: NED dates are in integer years, assume source imagery was flown in late summer for mountains
    if datefield == 'S_DATE_CLN':
        date_a += 0.75
    return date_a 
Example 6
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 5 votes vote down vote up
def ds2mask(v_ds, r_ds):
    lyr = v_ds.GetLayer()
    mask = lyr2mask(lyr, r_ds)
    lyr = None
    return mask

#Rasterize ogr layer to binary mask (core functionality for gdal.RasterizeLayer) 
Example 7
Project: firedpy   Author: earthlab   File: functions.py    MIT License 5 votes vote down vote up
def rasterize(src, dst, attribute, resolution, crs, extent, all_touch=False,
              na=-9999):

    # Open shapefile, retrieve the layer
    src_data = ogr.Open(src)
    layer = src_data.GetLayer()

    # Use transform to derive coordinates and dimensions
    xmin = extent[0]
    ymin = extent[1]
    xmax = extent[2]
    ymax = extent[3]

    # Create the target raster layer
    cols = int((xmax - xmin)/resolution)
    rows = int((ymax - ymin)/resolution) + 1
    trgt = gdal.GetDriverByName("GTiff").Create(dst, cols, rows, 1,
                                gdal.GDT_Float32)
    trgt.SetGeoTransform((xmin, resolution, 0, ymax, 0, -resolution))

    # Add crs
    refs = osr.SpatialReference()
    refs.ImportFromWkt(crs)
    trgt.SetProjection(refs.ExportToWkt())

    # Set no value
    band = trgt.GetRasterBand(1)
    band.SetNoDataValue(na)

    # Set options
    if all_touch is True:
        ops = ["-at", "ATTRIBUTE=" + attribute]
    else:
        ops = ["ATTRIBUTE=" + attribute]

    # Finally rasterize
    gdal.RasterizeLayer(trgt, [1], layer, options=ops)

    # Close target an source rasters
    del trgt
    del src_data 
Example 8
Project: surfclass   Author: Kortforsyningen   File: rasterio.py    MIT License 4 votes vote down vote up
def read_2d(self, geom):
        """Reads part of the raster into a 2D MaskedArray with a mask marking cells outside the polygon.

        This is suitable for doing analysis on cells inside a given geometry object.

        Note: The mask marks cells outside the geometry. This means that cells inside the geometry are NOT masked
            even if they are equal to the raster nodata value.

        Args:
            geom (osgeo.ogr.Geometry): OGR Geometry object

        Raises:
            TypeError: If geometry is not an `osgeo.ogr.Geometry`
            ValueError: If bbox of the geometry is entirely or partly outside raster coverage.

        Returns:
            [numpy.ma.maskedArray]: A masked array where cells outside the geometry are masked.

        """
        if not isinstance(geom, ogr.Geometry):
            raise TypeError("Must be OGR geometry")
        mem_type = geom.GetGeometryType()
        ogr_env = geom.GetEnvelope()
        geom_bbox = Bbox(ogr_env[0], ogr_env[2], ogr_env[1], ogr_env[3])
        window = self.bbox_to_pixel_window(geom_bbox)
        if window[2] <= 0 or window[3] <= 0:
            return np.ma.empty(shape=(0, 0))
        src_array = self.read_raster(window=window, masked=False)
        # calculate new geotransform of the feature subset
        new_gt = self.window_geotransform(window)

        # Create a temporary vector layer in memory
        mem_ds = self._ogr_mem_drv.CreateDataSource("out")
        mem_layer = mem_ds.CreateLayer("mem_lyr", geom.GetSpatialReference(), mem_type)
        mem_feature = ogr.Feature(mem_layer.GetLayerDefn())
        mem_feature.SetGeometry(geom)
        mem_layer.CreateFeature(mem_feature)

        # Rasterize the feature
        mem_raster_ds = self._gdal_mem_drv.Create(
            "", window[2], window[3], 1, gdal.GDT_Byte
        )
        mem_raster_ds.SetProjection(self._ds.GetProjection())
        mem_raster_ds.SetGeoTransform(new_gt)
        # Burn 1 inside our feature
        gdal.RasterizeLayer(mem_raster_ds, [1], mem_layer, burn_values=[1])
        rasterized_array = mem_raster_ds.ReadAsArray()

        # Mask the source data array with our current feature mask
        masked = np.ma.MaskedArray(src_array, mask=np.logical_not(rasterized_array))
        return masked 
Example 9
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 4 votes vote down vote up
def rasterize(self, src, dst, attribute, all_touch=False,
                  epsg=4326, na=-9999):
        '''
        It seems to be unreasonably involved to do this in Python compared to
        the command line.
        '''
        resolution = self.resolution

        # Open shapefile, retrieve the layer
        src_data = ogr.Open(src)
        layer = src_data.GetLayer()
        extent = layer.GetExtent()  # This is wrong

        # Create the target raster layer
        xmin, xmax, ymin, ymax = extent
        cols = int((xmax - xmin)/resolution)
        rows = int((ymax - ymin)/resolution)
        trgt = gdal.GetDriverByName('GTiff').Create(dst, cols, rows, 1,
                                   gdal.GDT_Float32)
        trgt.SetGeoTransform((xmin, resolution, 0, ymax, 0, -resolution))

        # Add crs
        refs = osr.SpatialReference()
        refs.ImportFromEPSG(epsg)
        trgt.SetProjection(refs.ExportToWkt())

        # Set no value
        band = trgt.GetRasterBand(1)
        band.SetNoDataValue(na)

        # Set options
        if all_touch is True:
            ops = ['-at', 'ATTRIBUTE=' + attribute]
        else:
            ops = ['ATTRIBUTE=' + attribute]

        # Finally rasterize
        gdal.RasterizeLayer(trgt, [1], layer, options=ops)

        # Close target an source rasters
        del trgt
        del src_data 
Example 10
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 4 votes vote down vote up
def rasterize(self, src, dst, attribute, all_touch=False,
                  epsg=4326, na=-9999):
        '''
        It seems to be unreasonably involved to do this in Python compared to
        the command line.
        ''' 
        resolution = self.resolution

        # Open shapefile, retrieve the layer
        src_data = ogr.Open(src)
        layer = src_data.GetLayer()
        extent = layer.GetExtent()  # This is wrong

        # Create the target raster layer
        xmin, xmax, ymin, ymax = extent
        cols = int((xmax - xmin)/resolution)
        rows = int((ymax - ymin)/resolution)
        trgt = gdal.GetDriverByName('GTiff').Create(dst, cols, rows, 1,
                                   gdal.GDT_Float32)
        trgt.SetGeoTransform((xmin, resolution, 0, ymax, 0, -resolution))
    
        # Add crs
        refs = osr.SpatialReference()
        refs.ImportFromEPSG(epsg)
        trgt.SetProjection(refs.ExportToWkt())
    
        # Set no value
        band = trgt.GetRasterBand(1)
        band.SetNoDataValue(na)
    
        # Set options
        if all_touch is True:
            ops = ['-at', 'ATTRIBUTE=' + attribute]
        else:
            ops = ['ATTRIBUTE=' + attribute]

        # Finally rasterize
        gdal.RasterizeLayer(trgt, [1], layer, options=ops)

        # Close target raster
        trgt = None
        src_data = None 
Example 11
Project: dask-geomodeling   Author: nens   File: utils.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def create_dataset(array, geo_transform=None, projection=None, no_data_value=None):
    """
    Create and return a gdal dataset.

    :param array: A numpy array.
    :param geo_transform: 6-tuple of floats
    :param projection: wkt projection string
    :param no_data_value: integer or float

    This is the fastest way to get a gdal dataset from a numpy array, but
    keep a reference to the array around, or a segfault will occur. Also,
    don't forget to call FlushCache() on the dataset after any operation
    that affects the array.
    """
    # prepare dataset name pointing to array
    datapointer = array.ctypes.data
    bands, lines, pixels = array.shape
    datatypecode = gdal_array.NumericTypeCodeToGDALTypeCode(array.dtype.type)
    datatype = gdal.GetDataTypeName(datatypecode)
    bandoffset, lineoffset, pixeloffset = array.strides
    # if projection is wrong there will be a segfault in RasterizeLayer
    projection = osr.GetUserInputAsWKT(str(projection))

    dataset_name_template = (
        "MEM:::"
        "DATAPOINTER={datapointer},"
        "PIXELS={pixels},"
        "LINES={lines},"
        "BANDS={bands},"
        "DATATYPE={datatype},"
        "PIXELOFFSET={pixeloffset},"
        "LINEOFFSET={lineoffset},"
        "BANDOFFSET={bandoffset}"
    )
    dataset_name = dataset_name_template.format(
        datapointer=datapointer,
        pixels=pixels,
        lines=lines,
        bands=bands,
        datatype=datatype,
        pixeloffset=pixeloffset,
        lineoffset=lineoffset,
        bandoffset=bandoffset,
    )

    # access the array memory as gdal dataset
    dataset = gdal.Open(dataset_name, gdal.GA_Update)

    # set additional properties from kwargs
    if geo_transform is not None:
        dataset.SetGeoTransform(geo_transform)
    if projection is not None:
        dataset.SetProjection(projection)
    if no_data_value is not None:
        for i in range(len(array)):
            dataset.GetRasterBand(i + 1).SetNoDataValue(no_data_value)

    return dataset 
Example 12
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 4 votes vote down vote up
def shp2array(shp_fn, r_ds=None, res=None, extent=None, t_srs=None):
    """Rasterize input shapefile to match existing raster Dataset (or specified res/extent/t_srs)
    """
    shp_ds = ogr.Open(shp_fn)
    lyr = shp_ds.GetLayer()
    #This returns xmin, ymin, xmax, ymax
    shp_extent = lyr_extent(lyr)
    shp_srs = lyr.GetSpatialRef()
    # dst_dt = gdal.GDT_Byte
    ndv = 0
    if r_ds is not None:
        r_extent = ds_extent(r_ds)
        res = get_res(r_ds, square=True)[0] 
        if extent is None:
            extent = r_extent
        r_srs = get_ds_srs(r_ds)
        r_geom = ds_geom(r_ds)
        # dst_ns = r_ds.RasterXSize
        # dst_nl = r_ds.RasterYSize
        #Convert raster extent to shp_srs
        cT = osr.CoordinateTransformation(r_srs, shp_srs)
        r_geom_reproj = geom_dup(r_geom)
        r_geom_reproj.Transform(cT)
        r_geom_reproj.AssignSpatialReference(t_srs)
        lyr.SetSpatialFilter(r_geom_reproj)
        #lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))
    else:
        #TODO: clean this up
        if res is None:
            sys.exit("Must specify input res")
        if extent is None:
            print("Using input shp extent")
            extent = shp_extent
    if t_srs is None:
        t_srs = r_srs
    if not shp_srs.IsSame(t_srs):
        print("Input shp srs: %s" % shp_srs.ExportToProj4())
        print("Specified output srs: %s" % t_srs.ExportToProj4())
        out_ds = lyr_proj(lyr, t_srs)
        outlyr = out_ds.GetLayer()
    else:
        outlyr = lyr
    #outlyr.SetSpatialFilter(r_geom)
    m_ds = mem_ds(res, extent, srs=t_srs, dtype=gdal.GDT_Byte)
    b = m_ds.GetRasterBand(1)
    b.SetNoDataValue(ndv)
    gdal.RasterizeLayer(m_ds, [1], outlyr, burn_values=[1])
    a = b.ReadAsArray()
    a = ~(a.astype('Bool'))
    return a 
Example 13
Project: wa   Author: wateraccounting   File: raster_conversions.py    Apache License 2.0 4 votes vote down vote up
def Vector_to_Raster(Dir, shapefile_name, reference_raster_data_name):
    """
    This function creates a raster of a shp file

    Keyword arguments:
    Dir --
        str: path to the basin folder
    shapefile_name -- 'C:/....../.shp'
        str: Path from the shape file
    reference_raster_data_name -- 'C:/....../.tif'
        str: Path to an example tiff file (all arrays will be reprojected to this example)
    """

    from osgeo import gdal, ogr

    geo, proj, size_X, size_Y=Open_array_info(reference_raster_data_name)

    x_min = geo[0]
    x_max = geo[0] + size_X * geo[1]
    y_min = geo[3] + size_Y * geo[5]
    y_max = geo[3]
    pixel_size = geo[1]

    # Filename of the raster Tiff that will be created
    Dir_Basin_Shape = os.path.join(Dir,'Basin')
    if not os.path.exists(Dir_Basin_Shape):
        os.mkdir(Dir_Basin_Shape)

    Basename = os.path.basename(shapefile_name)
    Dir_Raster_end = os.path.join(Dir_Basin_Shape, os.path.splitext(Basename)[0]+'.tif')

    # Open the data source and read in the extent
    source_ds = ogr.Open(shapefile_name)
    source_layer = source_ds.GetLayer()

    # Create the destination data source
    x_res = int(round((x_max - x_min) / pixel_size))
    y_res = int(round((y_max - y_min) / pixel_size))

    # Create tiff file
    target_ds = gdal.GetDriverByName('GTiff').Create(Dir_Raster_end, x_res, y_res, 1, gdal.GDT_Float32, ['COMPRESS=LZW'])
    target_ds.SetGeoTransform(geo)
    srse = osr.SpatialReference()
    srse.SetWellKnownGeogCS(proj)
    target_ds.SetProjection(srse.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    target_ds.GetRasterBand(1).SetNoDataValue(-9999)
    band.Fill(-9999)

    # Rasterize the shape and save it as band in tiff file
    gdal.RasterizeLayer(target_ds, [1], source_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
    target_ds = None

    # Open array
    Raster_Basin = Open_tiff_array(Dir_Raster_end)

    return(Raster_Basin) 
Example 14
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    GNU General Public License v3.0 4 votes vote down vote up
def clipRaster(raster, newRaster, vector):
    raster = gdal.Open(raster)
    
    vect = ogr.Open(vector)
    lyr = vect.GetLayer()
    ext = lyr.GetExtent()
    
    gTrans = raster.GetGeoTransform()
    #get the x start of the left most pixel
    xlP = int((ext[0] - gTrans[0])/gTrans[1])*gTrans[1] - abs(gTrans[0])
    #get the x end of the right most pixel
    xrP = math.ceil((ext[1] - gTrans[0])/gTrans[1])*gTrans[1] - abs(gTrans[0])
    #get the y start of the upper most pixel
    yuP = abs(gTrans[3]) - int((gTrans[3] - ext[3])/gTrans[5])*gTrans[5]
    #get the y end of the lower most pixel
    ylP = abs(gTrans[3]) - math.floor((gTrans[3] - ext[2])/gTrans[5])*gTrans[5]
        
    gdal.Translate('tmp.tif', raster, projWin = [xlP, yuP, xrP, ylP])
    del raster
    tRas = gdal.Open('tmp.tif')
    band = tRas.GetRasterBand(1)
    noDat = band.GetNoDataValue()
    # store a copy before rasterize
    fullRas = band.ReadAsArray().astype(numpy.float)
    
    gdal.RasterizeLayer(tRas, [1], lyr, None, None, [-9999], ['ALL_TOUCHED=TRUE']) # now tRas is modified
    
    finRas = tRas.GetRasterBand(1).ReadAsArray().astype(numpy.float)

    for i, row in enumerate(finRas):
        for j, col in enumerate(row):
            if col == -9999.:
                finRas[i, j] = fullRas[i, j]
            else:
                finRas[i, j] = noDat

    array2raster(newRaster, 'tmp.tif', finRas, noDat, "float32")
    os.remove('tmp.tif')
    del fullRas, finRas, band, tRas

# create a reference raster with random values    
# create a reference raster with random values 
Example 15
Project: chorospy   Author: spyrostheodoridis   File: bioFunc.py    GNU General Public License v3.0 4 votes vote down vote up
def makeDensityRaster(speciesOcc, inVector, pixelSize, outRas, noData):
    srcVector = ogr.Open(inVector)
    srcLayer = srcVector.GetLayer()
    srs = srcLayer.GetSpatialRef()
    # if the layer is not wgs84
    if srs.GetAttrValue("AUTHORITY", 1) != '4326':
        print('Layer projection should be WGS84!')
        return

    xMin, xMax, yMin, yMax = srcLayer.GetExtent()

    # Create the destination data source
    xRes = int((xMax - xMin) / pixelSize)
    yRes = int((yMax - yMin) / pixelSize)
    targetRas = gdal.GetDriverByName('GTiff').Create(outRas, xRes, yRes, 1, 6) # 6 == float
    targetRas.SetGeoTransform((xMin, pixelSize, 0, yMax, 0, -pixelSize))
    band = targetRas.GetRasterBand(1)
    band.SetNoDataValue(noData)

    # Rasterize clips the raster band
    gdal.RasterizeLayer(targetRas, [1], srcLayer, None, None, [0], ['ALL_TOUCHED=TRUE'])
    #os.remove(outRas)
    g = band.ReadAsArray()

    for point in speciesOcc.iterrows():
        xi = int((point[1]['x'] - xMin) / pixelSize)
        yi = int((point[1]['y'] - yMax) / -pixelSize)

        try:
            if g[yi,xi] != noData:
                g[yi,xi] += 1
            else:
                print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
        except:
            print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
            pass


    band.WriteArray(g)
    targetRasSRS = osr.SpatialReference()
    targetRasSRS.ImportFromWkt(srs.ExportToWkt())
    targetRas.SetProjection(targetRasSRS.ExportToWkt())
    band.FlushCache() 
Example 16
Project: watools   Author: wateraccounting   File: raster_conversions.py    Apache License 2.0 3 votes vote down vote up
def Vector_to_Raster(Dir, shapefile_name, reference_raster_data_name):
    """
    This function creates a raster of a shp file

    Keyword arguments:
    Dir --
        str: path to the basin folder
    shapefile_name -- 'C:/....../.shp'
        str: Path from the shape file
    reference_raster_data_name -- 'C:/....../.tif'
        str: Path to an example tiff file (all arrays will be reprojected to this example)
    """

    from osgeo import gdal, ogr

    geo, proj, size_X, size_Y=Open_array_info(reference_raster_data_name)

    x_min = geo[0]
    x_max = geo[0] + size_X * geo[1]
    y_min = geo[3] + size_Y * geo[5]
    y_max = geo[3]
    pixel_size = geo[1]

    # Filename of the raster Tiff that will be created
    Dir_Basin_Shape = os.path.join(Dir,'Basin')
    if not os.path.exists(Dir_Basin_Shape):
        os.mkdir(Dir_Basin_Shape)

    Basename = os.path.basename(shapefile_name)
    Dir_Raster_end = os.path.join(Dir_Basin_Shape, os.path.splitext(Basename)[0]+'.tif')

    # Open the data source and read in the extent
    source_ds = ogr.Open(shapefile_name)
    source_layer = source_ds.GetLayer()

    # Create the destination data source
    x_res = int(round((x_max - x_min) / pixel_size))
    y_res = int(round((y_max - y_min) / pixel_size))

    # Create tiff file
    target_ds = gdal.GetDriverByName('GTiff').Create(Dir_Raster_end, x_res, y_res, 1, gdal.GDT_Float32, ['COMPRESS=LZW'])
    target_ds.SetGeoTransform(geo)
    srse = osr.SpatialReference()
    srse.SetWellKnownGeogCS(proj)
    target_ds.SetProjection(srse.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    target_ds.GetRasterBand(1).SetNoDataValue(-9999)
    band.Fill(-9999)

    # Rasterize the shape and save it as band in tiff file
    gdal.RasterizeLayer(target_ds, [1], source_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
    target_ds = None

    # Open array
    Raster_Basin = Open_tiff_array(Dir_Raster_end)

    return(Raster_Basin)