4

I have a shapefile showing the locations of some cities and I am trying to extract pixel values from a RGB image at the location of points and store it in new columns in the point shapefile. I can print the extracted pixels values, but once I want to store them in new columns in the shapefile, I receive an error. can anyone help me to solve this using rasterio and geopandas?

Here is the code:

import geopandas as gpd
import rasterio as rio
nbands=3
cities=gpd.read_file('path to the point shapefile/cities.shp')
fields=['band1','band2','band3']
cities[fields]=0

for index, row in cities.iterrows():
    cities = row['cities']
    x= row['geometry'].x
    y= row['geometry'].y
    
    raster = rio.open('path to RGB image/stacked_RGB.tif') 
    raster_data= raster.read()
    for i in range(nbands):
        band = raster_data[:,:,i]
        rowIndex,colIndex = raster.index(x,y)
        print(cities + ': '+ str(raster_data[i][rowIndex,colIndex])) #up to here the code works fine
        cities[['band1','band2','band3']].loc[index]=raster_data[i][rowIndex,colIndex]

error:

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  field_samples[['band1','band2','band3']].iloc[index]=raster_data[i][rowIndex,colIndex]
/run/nvme/job_12710378/tmp/ipykernel_144365/871553810.py:14: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame
1

1 Answer 1

4

You're overcomplicating things. Just use the rasterio dataset.sample method.

E.g

import geopandas as gpd
import rasterio as rio

cities = gpd.read_file('cities.shp')
stacked_RGB = 'stacked_RGB.tif'

with rio.open(stacked_RGB) as raster: # Only open raster once, not every loop iteration
                                      # and definitely don't read() the raster every iteration
    fields = [f'band{i+1}' for i in range(raster.count)]  # calculate output fields instead of hardcoding
    cities[fields] = list(raster.sample(zip(cities['geometry'].x, cities['geometry'].y)))
    print(cities.head())

Here's an alternate implementation from the geopandas docs.

1
  • Great answer!it works fine.
    – rayan
    Aug 12, 2022 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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