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, geoT+nx*geoT) max_x = max(geoT, geoT+nx*geoT) min_y = min(geoT, geoT + geoT[-1]*ny) max_y = max(geoT, geoT + 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, 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