2

I want average tree cover for every census block. I have a file of census blocks, and I have a USFS tree cover file from USFS Tree Canopy Cover Datasets that looks like this:

over

I'm trying to use zonal statistics but the files are too big- even though my computer is very powerful QGIS only uses <3GB of RAM for some reason.

Is there a way to speed this up?

I tried using R's exact_extractr but I got all NaNs.

rasta = raster('canopy.tif')
poly3 = st_read("cbg.shp")
foliage <- exact_extract(rasta, poly3, 'mean')
4
  • In what way are the files "too big"? Do you get errors loading them into R? What's the resolution and size? What's the size of the census blocks compared to the raster cell size? How many census blocks are there? Can you run this for one census block? Can you subdivide into states or counties to make it easier? – Spacedman Sep 1 '20 at 11:18
  • 1
    I can't see why you'd get all NaNs with exactextractr (I'm the author). Can you share the specific inputs you're using? – dbaston Sep 2 '20 at 0:55
  • one of them is the first CONUS zip link at the USFS site I linked to above. The other is a gdb arcgis.com/home/item.html?id=1c924a53319a491ab43d5cb1d55d8561 – Ben Hendel Sep 2 '20 at 3:36
  • You're seeing NaNs because the first block groups in the file are from Hawaii, which is not covered by the CONUS raster. The results aren't "all NaN", but R only displays the first 1000 results by default, and there are 219,000 block groups. (3,940 of them are NaN). – dbaston Sep 3 '20 at 1:34
2

You could try using the GRASS GIS tool, v.rast.stats to see it is more efficient. In order to get set-up in GRASS with your Location and Mapset, please review the tutorial here and associated 'Getting Started docs'. Once inside your GRASS location,you could run the following commands:

##import your polygon
v.in.ogr input=cbg.shp output=cbg

##import your raster
r.in.gdal input=canopy.tif output=canopy

##calculate univariate statistics for each polygon, including mean
v.rast.stats map=cbg raster=canopy column_prefix=stat_ method=mean

You can also leave the method empty to calculate a range of default statistics.

Alternatively, you could try to run the tool from within QGIS, see screenshot below:

enter image description here Regarding your second question concerning extract_exact, perhaps you need to remove NAs, have you tried:

rasta = raster('canopy.tif')
poly3 = st_read("cbg.shp")
foliage <- exact_extract(rasta, poly3, 'mean', na.rm=TRUE)

More details here

While the approaches outlined above are valid, you may still need to consider tiling (e.g. gdal_retile.py) your dataset in order to reduce the associated file size of your raster.

4
  • what is a GRASS location? Is it in QGIS? R? – Ben Hendel Sep 1 '20 at 10:18
  • a GRASS location is specific to GRASS-GIS (see edited answer above). Alternatively run the tool in QGIS if it is easier. – dmci Sep 1 '20 at 10:26
  • I tried na.rm = true and I got Error in .local(x, y, ...) : exact_extract was called with a named summary operation that does not accept additional arguments ... In addition: Warning message: In .local(x, y, ...) : Polygons transformed to raster CRS (EPSG:NA) – Ben Hendel Sep 1 '20 at 21:28
  • And when I use v.rast.stats I get all NA values for the mean... – Ben Hendel Sep 1 '20 at 21:40
0

I ended up using regular Zonal Statistics in QGIS, I let it run for four hours but it worked. Other methods just gave me NA as the mean

1
  • If you are running QGIS > 3.14, I'd recommend to have a look at the command line tool qgis_process. Thus you can use qgis tools from scripts, e.g. from R. – loki Jan 31 at 14:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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