"Spatial join" between raster and shapefille- give raster the shapefile values with python

Question

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):

the shapefile has column which contain data about each polygon:

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

Answer 1

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.