Assigning an ID to each pixel which overlaps with polygon shapefile (Python)

Question

I calculated NDVI values for the last 24 years and created Pandas Dataframe which every cell includes columns with XY coordinates and 24 NDVI values.

I have several polygon shapefiles which overlap some of the pixels.

I would like to find pixels which overlap with shapefile polygon and assign some ID (for example 1 if overlapped and 0 if not overlapped) in a separate column in geodataframe.

I found some ArcMap tools which can be used here: Using polygons to apply unique value to all raster cells within those polygons? and here: Assign raster cells presence/absence values depending on whether polygon overlaps with cells.

Is there any python code which can do this?

gdfdata.head()
out[]
id  ndvi2017    geometry    ndvi1994    ndvi1995    ndvi1996    ndvi1997    ndvi1998    ndvi1999    ndvi2000    ... ndvi2007    ndvi2008    ndvi2009    ndvi2010    ndvi2011    ndvi2012    ndvi2013    ndvi2014    ndvi2015    ndvi2016
0   1   0.377508    POINT (1553412.5 -3946412.5)    0.440593    0.501229    0.614488    0.457669    0.516372    0.652276    0.574935    ... 0.442016    0.429176    0.253007    0.235713    0.259660    0.334798    0.301211    0.390102    0.355966    0.498201
1   2   0.450151    POINT (1553437.5 -3946412.5)    0.438102    0.508961    0.608286    0.462507    0.520403    0.690018    0.590940    ... 0.456732    0.469569    0.332611    0.325883    0.373628    0.476741    0.378534    0.552504    0.437974    0.615305
2   3   0.412462    POINT (1553437.5 -3946437.5)    0.426529    0.500758    0.589582    0.453873    0.513041    0.663890    0.563998    ... 0.441368    0.464280    0.351082    0.310373    0.395219    0.466424    0.383781    0.504997    0.440410    0.620843
3   4   0.454601    POINT (1553462.5 -3946387.5)    0.454064    0.513976    0.610625    0.475956    0.538307    0.711814    0.604684    ... 0.452001    0.457405    0.364040    0.442932    0.421774    0.524144    0.418083    0.644018    0.464420    0.662594
4   5   0.489210    POINT (1553462.5 -3946412.5)    0.452784    0.531119    0.610266    0.488797    0.543849    0.710351    0.601440    ... 0.468247    0.512874    0.443819    0.471668    0.493467    0.576141    0.488978    0.677011    0.493551    0.690452

shapefile:

shp_fence = gpd.read_file('fence.shp')
print(shp_fence)  

out[]
  Id                                           geometry
0   1  POLYGON ((1555435.047127967 -3944370.721410914...
1   2  POLYGON ((1554617.469131603 -3944278.199067593...
2   3  POLYGON ((1556472.315180207 -3945250.625474782...
3   4  POLYGON ((1555910.622499535 -3945103.390798647...
4   5  POLYGON ((1555526.40089912 -3946411.56900058, ...
5   6  POLYGON ((1556258.591632321 -3946544.835142805...
6   7  POLYGON ((1555270.818279916 -3949387.08823332,...
7   8  POLYGON ((1556123.778012351 -3949736.517403783...

the same coordinate system:

print(gdfdata.crs)
print(shp_fence.crs)
out[]
{'init': 'epsg:3577'}
{'init': 'epsg:3577'}

print(type(gdfdata))
print(type(shp_fence))
out[]
<class 'geopandas.geoseries.GeoSeries'>
<class 'geopandas.geodataframe.GeoDataFrame'>

They overlap (red color are polygons and blue are points):

I used GeoSeries.within() to find points within my shapefile:

gdfdata['fenced'] = gdfdata.within(shp_fence)

UPDATE: Warning:

 C:\Users\u6262380\AppData\Local\Continuum\anaconda3\lib\site-packages\ipykernel_launcher.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

 See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy

However, this code runs and returns a new column with FALSE values only although I have a lot of points inside polygon which should return TRUE values. I wonder what the 'warning' above mean, I suppose that warning is causing a problem.

UPDATE:

It is suggested to iterate and check each point if it is in polygon here https://automating-gis-processes.github.io/2016/Lesson3-point-in-polygon.html?highlight=within I tried to do that with following code:

 gdfdata['fence'] = ''
 for row in gdfdata.geometry:
     gdfdata['fence'] = row.within(shp_fence)

But getting following error:

AttributeError: 'GeoDataFrame' object has no attribute '_geom'

Can someone tell me what is wrong?

Answer 1

I found a solution here. Basically it creates new column for each polygon and checks if the point is inside that polygon. If the point is inside polygon, the point is assigned True value, if outside False value is given. I found it the most elegant way of finding the points and corresponding polygons.

shp_fence=shp_fence.geometry
pnts = gdfdata.assign(**{key: gdfdata.within(geom) for key, geom in shp_fence.items()})

print(pnts)