Calculate solar irradiance inside a geometry trough rasterstats

Question

I am using the rasterstats python module zonal_stats to calculate the mean solar irradiance value (taken from a raster) inside each land particle (contained in a shapefile).

If I run print(irradiance.read().shape) I get (1, 2852, 2425).

Then, I do:

with rioxarray.open_rasterio("C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif") as src:
    stats = zonal_stats(particles,src)

And I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_14588\3822157982.py in <module>
      1 with rioxarray.open_rasterio("C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif") as src:
----> 2     stats = zonal_stats(particles,src)

~\anaconda3\lib\site-packages\rasterstats\main.py in zonal_stats(*args, **kwargs)
     38     The only difference is that ``zonal_stats`` will
     39     return a list rather than a generator."""
---> 40     return list(gen_zonal_stats(*args, **kwargs))
     41 
     42 

~\anaconda3\lib\site-packages\rasterstats\main.py in gen_zonal_stats(vectors, raster, layer, band, nodata, affine, stats, all_touched, categorical, category_map, add_stats, zone_func, raster_out, prefix, geojson_out, boundless, **kwargs)
    163         band = band_num
    164 
--> 165     with Raster(raster, affine, nodata, band) as rast:
    166         features_iter = read_features(vectors, layer)
    167         for _, feat in enumerate(features_iter):

~\anaconda3\lib\site-packages\rasterstats\io.py in __init__(self, raster, affine, nodata, band)
    272             self.nodata = nodata
    273         else:
--> 274             self.src = rasterio.open(raster, "r")
    275             self.affine = guard_transform(self.src.transform)
    276             self.shape = (self.src.height, self.src.width)

~\anaconda3\lib\site-packages\rasterio\env.py in wrapper(*args, **kwds)
    449 
    450         with env_ctor(session=session):
--> 451             return f(*args, **kwds)
    452 
    453     return wrapper

~\anaconda3\lib\site-packages\rasterio\__init__.py in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)
    244             or isinstance(fp, (os.PathLike, MemoryFile, FilePath))
    245         ):
--> 246             raise TypeError("invalid path or file: {0!r}".format(fp))
    247     if mode and not isinstance(mode, str):
    248         raise TypeError("invalid mode: {0!r}".format(mode))

TypeError: invalid path or file: <xarray.DataArray (band: 1, y: 3000, x: 2000)>
[6000000 values with dtype=int32]
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 7.6e+05 7.6e+05 7.6e+05 ... 7.8e+05 7.8e+05 7.8e+05
  * y            (y) float64 4.09e+06 4.09e+06 4.09e+06 ... 4.06e+06 4.06e+06
    spatial_ref  int32 0
Attributes: (12/15)
    AREA_OR_POINT:           Area
    BandName:                Band_1
    RepresentationType:      ATHEMATIC
    STATISTICS_COVARIANCES:  8498387208.982427
    STATISTICS_MAXIMUM:      1552768
    STATISTICS_MEAN:         1286479.1959047
    ...                      ...
    STATISTICS_SKIPFACTORY:  1
    STATISTICS_STDDEV:       92186.69757065
    _FillValue:              -1
    scale_factor:            1.0
    add_offset:              0.0
    long_name:               Band_1

It seems an error related to the use of float64 or int32, but I can't understand how to make it work. Is there a way to do it? And, in general, what is the easiest way to calculate the average value of some pixels inside polygons? (Given in the same shapefile I have many polygons on which I have to perform the mean of the pixels)

Answer 1

The second parameter of the zonal_stats function can just be the file path to the raster file, so something like this should work:

raster_path = "C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif"
stats = zonal_stats(particles, raster_path)