2

I am trying to use a for loop which would produce the zonal statistics of multiple rasters at shapefile points and then add each raster stat to the shapefile.

I've gotten it to add a column, but I don't know how to make the column change the name from 'max' to what the raster is or an abbreviation of it so it doesn't add the rest of the zonal_stats output.

import rasterstats as rs
from raster2xyz.raster2xyz import Raster2xyz
import os, fnmatch

inputFolder = "C:/585L/PROJECT/precip/clipped_2016/"
fire_ign = gpd.read_file("C:/585L/PROJECT/fire_points/fire_pts_clip2016_2021.shp")
fire_ign_slim = fire_ign[['FIRE_YEAR', 'MTBS_FIRE_', 'DISCOVERY_', 'FIRE_SIZE', 'LATITUDE', 'LONGITUDE', 'FIPS_NAME', 'geometry']]

def findRasters (path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join (root, file)

            
for raster in findRasters (inputFolder, '*.bil'):
    stats = rs.zonal_stats(fire_ign_slim,
                      raster, 
                      stats = "max")
    statsdf = pd.DataFrame(stats)
    df_redux = pd.concat([statsdf, fire_ign_slim], axis = 1)

Results in:

max     FIRE_YEAR   MTBS_FIRE_  DISCOVERY_  FIRE_SIZE   LATITUDE    LONGITUDE   FIPS_NAME   geometry
616.276001  2016    CLARK   8/4/2016    2819.0  37.798889   -118.922500     Mono County     POINT (-118.92249 37.79888)

The clipped_2016 folder has 13, one for each month of the year and an annual, rasters in it though, and I would like to see the max of all 13 rasters with an appropriate name.

So it might read something like:

FIRE_YEAR   MTBS_FIRE_  DISCOVERY_  FIRE_SIZE   LATITUDE    LONGITUDE   FIPS_NAME   geometry    ppt_201601  ppt_201602  ppt_201603  ppt_201604
2016    CLARK   8/4/2016    2819.0  37.798889   -118.922500     Mono County     POINT (-118.92249 37.79888)     some number from raster     some number from raster     some number from raster     some number from raster

But I am stuck on figuring out how to take each zonal stat from each point from each raster and putting those numbers into a df with unique column named from the raster it was taken from.

And this is how my clipped_2016 folder looks if it is helpful: enter image description here

1 Answer 1

3

You can use os.path (basename and splitext) to get the raster name without path and extension and just add the max column to the dataframe with a new name: e.g.

import os.path

etc...

for raster in findRasters (inputFolder, '*.bil'):
    stats = rs.zonal_stats(fire_ign_slim,
                      raster, 
                      stats = "max")
    raster_name = os.path.splitext(os.path.basename(raster))[0]
    statsdf = pd.DataFrame(stats)
    fire_ign_slim[f"max_{raster_name}"] = statsdf["max"]

Or if you prefer concat, you can rename the max column:

for raster in findRasters (inputFolder, '*.bil'):
    stats = rs.zonal_stats(fire_ign_slim,
                      raster, 
                      stats = "max")
    raster_name = os.path.splitext(os.path.basename(raster))[0]
    statsdf = pd.DataFrame(data={'max': [1, 2, 3]}).rename(columns={"max": f"max_{raster_name}"})
    fire_ign_slim = pd.concat([fire_ign_slim, statsdf], axis=1)

Note that fire_ign_slim is the output from concat, not df_redux otherwise you are continually adding only one column.

3
  • Thank you, this works for renaming the column, but it still only adds the first raster in the folder to the dataframe when using 'df_redux = pd.concat([statsdf, fire_ign_slim], axis = 1)'. I was hoping each raster's zonal stats would be added to the dataframe but maybe I am looping the concat wrong, I don't know
    – sarakota
    Mar 7, 2023 at 21:15
  • 1
    @sarakota you only get one column added as you concat statsdf to fire_ign_slim and output to df_redux instead of fire_ign_slim
    – user2856
    Mar 7, 2023 at 21:57
  • Thank you very much, this solved my question, and thank you for explaining what I did wrong. That makes it easier for me to understand
    – sarakota
    Mar 7, 2023 at 22:02

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.