For anyone else that stumbles upon this you can do:
To get just one stat:
gdf['mean'] = pd.DataFrame(
To get all computed stats:
gdf = gdf.join(
You can use Zonal Raster Statistics (also called Zonal Grid) tool to count the number of pixel in each class. The tool will create a new CSV file with the number of pixels for each category.
The only required parameter is the classified image (Zonal Grid) that contains the categories you want to calculate.
In the following example, I have a land use ...
You can do this in R with the help of raster and rgdal library
Load the shapefile and raster
shp <- readOGR("shape file directory", "shapegile_name")
in.ras <- raster("raster file")
Extract the raster value over the shapefile
val <- extract(in.ras, shp)
Create a data frame to store the extracted value and the ...
You need to use layer.source() to get the paths of the rasters which is required by the zonal_stats module. Your code should look like:
from rasterstats import zonal_stats
layers = QgsMapLayerRegistry.instance().mapLayers().values()
for layer in layers:
stats = zonal_stats("/home/myshape.shp", layer.source())
You could run a benchmark test between QGIS and pktools and specifically pkextract (extract pixel values from raster image from a (vector or raster) sample) to see if it is faster when run directly from the command-line. The usage of the tool is outlined on the above-mentioned link, but in general terms it is:
pkextract -f 'ESRI Shapefile' -s vector_aoi....
The results of dataset.transform and dataset.affine are equals (see Migrating to Rasterio 1.0: affine.Affine() vs. GDAL-style geotransforms)
With rasterio version < 1 (affine = affine transform, transform = GDAL-style geotransforms,
with rasterio.open('slope.tif') as src:
affine = src.affine
The only way for knowing if both layers are aligned properly is putting them together at the Map Canvas. I did that (see next image), installed rasterstats python module with easy_install and, finally, ran a similar example as in your link.
>>>from rasterstats import zonal_stats
>>>stats = zonal_stats("/home/zeito/pyqgis_data/new_polygon....
Your raster DEM and vector road shapefile have different coordinate systems. They need a common coordinate reference system for the rasterstats package to compute the zonal statistics.
import geopandas as gpd
import rasterstats as rs
src = rasterio.open('example_dem.tif')
gdf = gpd.read_file('example_road.shp')
You might want to take a look at the SAGA Zonal Grid statistics. This is available in Processing if you have SAGA GIS installed. I've not tried this yet but it looks like it will do what you want.
According to the documentation this can create a contingency table (cross-tabulation) which includes the number of pixels ("n") in each category for each polygon. ...
Not a direct solution but a possible (and possibly tedious) workaround is to:
Use the Split vector layer tool on your "veg_class" column which creates a shapefile for each classification.
Then run the Zonal Statistics plugin (which allows you to choose which statistics you want calculated) or the tool from the Processing Toolbox or Raster menu on each ...
Notice this part of your code:
for raster in Lista_raster:
stats = zonal_stats("Pivo_4_wgs.shp", Lista_raster, # <---here
You're passing the entire list for each iteration of your loop. Try passing the the element of the iterator instead:
for raster in Lista_raster:
stats = zonal_stats("Pivo_4_wgs.shp", raster, ...
You're also only writing a ...
If you are doing linear modelling (you don't say what your model is, but this applies fairly widely) then no, you can't do this, it doesn't help.
You can't add collinear terms to a model. In other words, you can do this:
y ~ x1 + x2
y ~ x1 + x2 + 3*x1
and expect to get a term for x1 and 3*x1. (You can try this in R but you have to wrap the last ...
There are many raster calculators available. Depends, what software you are familiar with.
E.g. you can use GRASS GIS r.series module (the module is also accessible via QGIS toolbox).
As @mdsumner commented, you can also use R.
A programatic approach is to read data by Python Numpy and calculate desired statistics. Numpy methods are also implemented in ...
The function help of zonal_stats links to gen_zonal_stats since all the parameters you pass to raster_stats are passed to gen_zonal_stats.
I think the parameter you are looking for is nodata. The help for this parameter says: If raster is a GDAL source, this value overrides any NODATA value specified in the file’s metadata. If None, the file’s metadata’s ...
"If it is possible in Python, will weighted averaging that way be faster?"
See answer to next question. Collecting the aggregated stats along with the rasterizing under the hood is what takes the time. You can get your weighted average using the following answer. As far as what is faster, this I cannot comment on...if you're willing to re-code rasterstats,...