0

I have raster which has NDVI values .In addition, I have shapefile(geopandas) on top of this raster but it doesn't fully covers the raster (illustration):

enter image description here

the shapefile has column which contain data about each polygon:

enter image description here

I would like to do something like spatial join between them, to have this data (e.g X,Y,Z,A,C) stored for each pixel, so in the end i'lll have "new band" which has one band of NDVI and one with this letters data, so for example, pixel that doesn't fall inside the polygon gets value "null" or "nodata", and pixels that falls inside the polygon will get the value of the letter according to where they fall inside (if thery are inside polygon X they get value X, if they are inside polygon Z then they get Z ect.

yet I haven't found any way to do that, I though maybe to clip the raster with the shapefile and then give it the values but also then i'm not sure how to give it the different letters values.

My end goal: To have the shapefile data (which is stores in one column in the table) stored for each pixel, if the pixel is not inside the polygon - give no data value. Kind of spatial join for raster-pixel .

Edit: i'm using python with jupyter notebook .

I had the thought maybe to extract the raster to point and then try to do spatial join and back to raster but it's huge number of pixels

Edit2: I have tried Jose's solution to use this function:

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)

but when I try to apply this I get error:

    # Do the magic
    ds_dst= gdal.Rasterize(img, plots, attribute='Crop',
                        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

NameError: name 'srs' is not defined

1
  • @BERA Rasterio, geopandas and earthpy but open to new things
    – Reut
    Jun 18, 2020 at 9:43

1 Answer 1

2

You can do that quite easily directly with GDAL. The following function accomplishes what you want


import gdal # That's the only import you need, I think...

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

You need to call the above function with the raster and vector filenames as

rasterise_me("path/to/my/file.tif", "/path/to/my/shapefile.shp", "crop")

(need to change the filenames to reflect your data, and the attribute is set to crop, but you'd need to change it). Note that you cannot encode text, but integers, so if crop is a text variable, you'd need to translate it into a numbers and create a new field to rasterise on.

4
  • thank you for youranswer, I have tried it but when I trie the "do the magic" part I get error NameError: name 'srs' is not defined, I edit and added it in the original post. It's important to mention that the values are not numerical but it's crop's name, and the image has 4 bands.
    – Reut
    Jun 21, 2020 at 11:21
  • and tried it also with non string attribute but still the same error
    – Reut
    Jun 21, 2020 at 12:36
  • I think you need to copy and paste the whole function and then call it with your relevant filenames (I tested it with my data and it works). Also, if your shapefile field is text, you need to add a new field that encodes it as a number (e.g "wheat"->1, "barley"->2 etc)
    – Jose
    Jun 22, 2020 at 16:00
  • One thing concerning the error: check that your file has projection set up (e.g. check with gdalinfo on the console). If still not working, please maybe add a new focused question.
    – Jose
    Jan 17, 2022 at 14:51

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.