14

The below code worked for me QGis 1.8.0 You might modify this to accomodate multiple files with some loop.. from qgis.analysis import QgsZonalStatistics #specify polygon shapefile vector polygonLayer = QgsVectorLayer('F:/temp/zonalstat/zonePoly.shp', 'zonepolygons', "ogr") # specify raster filename rasterFilePath = 'F:/temp/zonalstat/raster1.tif' # ...


10

There is a bug that seems to correspond to what you're experiencing - it's registered as BUG-000084883 - The 'Ignore NoData in calculations' option in Zonal Statistics as Table tool {and Zonal Statistics tool} is not honored when checked off, producing incorrect results. It occurs with 10.3 and 10.2.2 but not 10.1. Did you try the tool with this version?


10

What you are looking for is COUNT, which is the frequency of the cells that you processed through Zonal Statistics. Sum, on the other hand, is the sum of cell values covered by your polygon. Overly simplistically, say, your cell values are 2,1,3,4,4 in this case COUNT is 5 and SUM is 14.


9

From Esri's support site : HowTo: Create points representing the highest or lowest elevations within polygon features Just replace the elevation raster by the Flow Accumulation raster. Identify the value of the highest elevation within each polygon feature using the Zonal Statistics tool: Open ArcMap and navigate to ArcToolbox > Spatial Analyst ...


9

It is a bug. Something terribly wrong with cell count. Correct mean (9.0452380952381) times correct number of non-empty cells (420) divided by 297 (that is a cell count reported by tool) results in 12.7912457912458. That is a wrong average reported by tool. Results of my own toy size grids test:


8

For anyone else that stumbles upon this you can do: To get just one stat: gdf['mean'] = pd.DataFrame( zonal_stats( vectors=gdf['geometry'], raster='raster.tif', stats='mean' ) )['mean'] To get all computed stats: gdf = gdf.join( pd.DataFrame( zonal_stats( vectors=gdf['geometry'], raster='...


7

With multiprocessing, for fastness! Has a little different output-formatting. #!/usr/bin/python import gdal, ogr, osr, numpy, sys from multiprocessing import Pool # Raster dataset input_value_raster = sys.argv[1] # Vector dataset(zones) input_zone_polygon = sys.argv[2] # Open data rast = gdal.Open(input_value_raster) shp = ogr.Open(input_zone_polygon) # ...


7

The Modifiable Aerial Unit Problem (MAUP) is a change of support issue associated with arbitrary aggregate units. Two classic examples are census tracks and wildlife game units. These have been found to be arbitrary political units and the underlying statistical response in demography acts independent of the unit. Because of this, the unit is not an accurate ...


7

QGIS first makes an initial pass, checking to see if the center of each raster cell is within the polygon. If fewer than two cell centers are within the polygon, it performs a vector-based intersection for all intersecting cells, whether their center is within the polygon or not, and computes a weight that is the fraction of each cell that is covered by the ...


7

If you wrote zonal in the Processing toolbox search, you will find many: The zonal Statistics is the second one from the top.


6

If you want to get zonal statistics for several features in one shapefile, you have to loop over the zonal_stats function. You can write the results of the loop for example to a dictionary. Below is the modified zonal_stats function together with a loop, looping over the input shapefile. As an output you get a Dictionary containing for each Feature ID the ...


6

there is another function called "zonal histogram" in ArcGIS (http://help.arcgis.com/en%20/arcgisdesktop/10.0/help/index.html#//009z000000w6000000.htm). With the table you can extract anything you want. It is however recommende to resample your nightime data before use in order to have an integer image. Don't forget to set the pixel size equal to the ...


6

The example data from Jeffrey library(raster) r <- raster(ncols=10, nrows=10) set.seed(0) x <- runif(ncell(r)) x[round(runif(25,1,100),digits=0)] <- NA r[] <- x cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)) cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)) polys <- SpatialPolygons(list(...


6

Use arcpy.env.overwriteOutput. arcpy.env.overwriteOutput = True # Execute ZonalStatisticsAsTable outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster, outTable, "NODATA", "MEAN")


6

I think you hit on your best option, which is to convert the raster to a vector and then intersect the result with your polygon layer. As a way of explanation regarding the frustration you're experiencing with Zonal Statistics (and actually this will also apply to your idea of cropping or "clipping" the raster), there is no alternative way for this to ...


6

QGIS and R's raster package use different methods to estimate zonal statistics. Briefly: QGIS compares the centroid of each raster cell to the polygon boundary, initially considering cells to be wholly within or outside of the polygon based on the centroid. However, if fewer than two cell centroids fall within the polygon, an exact vector-based calculation ...


5

For a count of "lit" pixels--or of any other kind of cell value for that matter--simply create a binary indicator grid. This is a grid with ones at the cells where the values are to be counted and zeros elsewhere: it should be clear that the zonal sum of these values counts the cells. To create an indicator, exploit the fact that true values will be stored ...


5

Before using Zonal Statistics make sure that you have same projections for both raster data and vector (shapefile) data. Sometimes different projection produce empty results. In the following example, I used WGS 84 for both raster data (SRTM Global) and a polygon vector data (Test.shp), and I didn't specify any column prefix: The final output is as follows:


5

When pixels are large comparing with polygons you better go "vector way". i.e. vectorize the raster tiles and then procede to a vector/vector intersection/computation. If you use the PostGIS Addons, you can do it like this: SELECT gt.id, (aws).geom, (aws).totalarea, (aws).weightedmean, FROM (SELECT id, ST_AreaWeightedSummaryStats(gv) ...


5

For calculating stats from raster cells intersected by line, you can use GRASS v.rast.stats directly from QGIS processing toolbox. It can calculate 13 different stats.


5

In R, you can extract the raster data for each polygon and then summarize it. First, lets create some data (FYI, you can read in a shapefile using raster::shapefile or rgdal::readOGR and a raster using raster::raster). library(raster) library(rgeos) r <- raster::raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, ymx=...


4

what I would recommend is to avoid zonal statistics when your zones are smaller than the pixel size, especially in your case where you have a regular grid. Instead, you should get the centroids of your polygons, then use the extract multiple value to point. There is an option for the interpolation.


4

This is explained in the Zonal Statistics help: If the zone input is a feature dataset with relatively small features, keep in mind that the resolution of the information needs to be appropriate relative to the resolution of the value raster. If the areas of single features are similar to or smaller than the area of single cells in the value raster, in the ...


4

I would use a zonal statistics analysis to investigate the relationship between crime rate and land price. Use Zonal Statistics as Table (Spatial Analyst) tool with these settings: in_zone_data = crime rate zone_field = label (field of crime rate) in_value_raster = land price statistics_type = all The result is a table. Use Excel to create a chart. ...


4

Begin with the elevation raster (the DSM). Compute focal minima and maxima using 2m radial neighborhoods, producing two new rasters. Subtract the original raster from these focal rasters. Take the absolute values of the differences. Compute the local maximum of the two absolute differences in (2): this gives the largest height deviation within 2 meters of ...


4

1) string formatting: It is a not a problem of QGIS or processing but a simple problem of string formatting in Python, look at Python String Format, for example. print "Var_{}".format(3) Var_3 # or old print "Var_%s" % (1) Var_1 Script example 1: for i in range(1,20): processing.runalg("qgis:zonalstatistics", "/folder/input/raster.tif",i,"/folder/...


4

Unfortunately, because ArcGIS's source code is not publicaly available, we cannot know for certain how ESRI treats boundary locations when you provide a vector zone layer input. However, as DanC points out above, it is very likely that there is some kind of internal vector-to-raster conversion that is taking place such that the vector zone layer is mapped ...


4

The Focal Statistics tool in Spatial Analyst can take care of the local maxima. Unfortunately, the size and shape of the neighborhood around the focal cell is fixed for the entire run. Thus, you will need to rerun Focal Statistics for each window size by adjusting the neighborhood parameter: Once you have one local maximum raster for each window size, you ...


4

We have one of these files on hand: ## dp is the root to our local data repository f <- file.path(dp, "data", "ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2.derived/pressure/air.mon.mean.nc") library(raster) b <- brick(f, level = 1) b On my system that uses the "ncdf4" package, but it could also use "ncdf". The global Pacific-view extent looks close, ...


4

With a multi-band raster input, ArcGIS will only process the first band: Multiband raster data When a multiband raster is used as input, most Spatial Analyst tools operate only on the first band. The exceptions are certain tools in the Multivariate and Extraction toolsets which do process each of the bands in a multiband input and can create a ...


Only top voted, non community-wiki answers of a minimum length are eligible