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