0

I have used rasterize function from GDAL in order to rasterize a shapefile that I have according to existing raster. The problem is that I'm not able to find or to show the resulting raster.

This is the rasterize function I am using:

def rasterise_me(raster, vector, attribute,
                fname_out="", format="MEM"):
    """Rasterises a vector dataset by attribute to match a given
    raster dataset. This functions allows for the raster and vector
    to have different projections, and will ensure that the output
    is consistent with the input raster.
    
    By default, it returns a handle to an open GDAL dataset that you
    can e.g. `ReadAsArray`. If you want to generate a  GTiff on disk,
    set format to `GTiff` and `fname_out` to a sensible filename.
    
    Parameters
    ----------
    raster: str
        The raster filaname used as input. It will not be overwritten.
    vector: str
        The vector filename
    attribute: str
        The attribute that you want to rasterize. Ideally, this is
        numeric.
    fname_out: str, optional
        The output filename.
    format: str, optional
        The output file format, such as GTiff, or whatever else GDAL
        understands
    """
    # Open input raster file. Need to do this to figure out
    # extent, projection & resolution.
    g = gdal.Open(raster) 
    geoT = g.GetGeoTransform()
    nx, ny = g.RasterXSize, g.RasterYSize 
    srs = g.GetProjection()
    min_x = min(geoT[0], geoT[0]+nx*geoT[1])
    max_x = max(geoT[0], geoT[0]+nx*geoT[1])
    min_y = min(geoT[3], geoT[3] + geoT[-1]*ny)
    max_y = max(geoT[3], geoT[3] + geoT[-1]*ny)
    # Reproject vector to match raster file
    vector_tmp = gdal.VectorTranslate("", vector, format="Memory",
                                    dstSRS=srs)
    # Do the magic
    ds_dst= gdal.Rasterize(fname_out, vector_tmp, attribute=attribute,
                        outputSRS=srs, xRes=geoT[1], yRes=geoT[-1],
                        outputBounds=[min_x, min_y, max_x, max_y],
                        format=format, outputType=gdal.GDT_Int32)
    return ds_dst

then I apply it on my data:

rasterize=rasterise_me(r"the/path/to/my/raster.tif", r"the/path/to/my/shape.shp", "class",fname_out='Rasterize.tif')

I couldn't find the results so I had tried to define the path also for fname_out like this:

rasterize=rasterise_me(r"the/path/to/my/raster.tif", r"the/path/to/my/shape.shp", "class",fname_out='this/is/the/path/for/my/future/raster/Rasterize.tif')

but still couldn't find it.

My endgoal is to be able to find my lost rasterize result raster.

0

1 Answer 1

1

You're writing a "MEM" raster, it's in MEMory. If you want to write to file, use a file based format like a GeoTIFF.

From the docstring in @Jose's rasterize_me function:

If you want to generate a GTiff on disk, set format to GTiff and fname_out to a sensible filename.

So do something like:

rasterize=rasterise_me(r"the/path/to/my/raster.tif", r"the/path/to/my/shape.shp", "class", fname_out='this/is/the/path/for/my/future/raster/Rasterize.tif', format="GTIFF")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.