1

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

enter image description here

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?

  • @BERA, ok, I will provide my code soon. – Sher Sep 6 '18 at 7:42
  • @ahmadhanb, I added the code, please activate my question. – Sher Sep 10 '18 at 5:18
  • 1
    @BERA, I added the code, please activate my question. – Sher Sep 10 '18 at 5:22
  • Can you provide the full minimal code to reproduce such error message? – Andre Silva Sep 10 '18 at 13:26
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.