1

I am trying to re-create my ArcGIS Zonal Stat model with R All my files will be 30m pixels; 18000 x 18000. I have used R two different ways to create the output and both create NoData values in the bottom portion of the output raster. My Arc GIS Zonal Stat Max output has values of zero & one. My R zonal Stat Max is mostly zero & one except for a portion of the image; ~~ 86 pixels from the bottom up, across the entire bottom, primarily come out as zero & No Data. I have included a picture of my output rasters showing the one section within those 86 pixels that actually write values of zero. I cannot just flip my NoData values to zero as there are other sections that need a value of one that are also getting set to NoData.

Python 2.7 with rpy2 created the clump file. To create the Zonal File, I needed to write an R script file then CALL R within python. I used the link below to write my R script.

http://www.guru-gis.net/efficient-zonal-statistics-using-r-and-gdal/

system.time(q <- reclassify(z, as.matrix(Zstat)[,1:2]))

Where z is my zonal file & Zstat is my zonal data frame.

After running line by line in R, I figured out the line of code above writing of the output file was taking the longest. 12 hours to complete or create the (q variable) on a 18000 x 18000 30m file PLUS NoData values in the bottom ~ 86 pixels of the output. I then changed the code to create a matrix from the data frame with code

I found here: https://stackoverflow.com/questions/52073767/joining-tables-to-a-raster-attribute-table-to-create-new-raster-with-the-reclass

system.time(reclass <- matrix(Zstat[,"vals"][as.vector(z)+1]))

This now takes 20 min & another 20 min to create (q - below)

system.time(q <- raster(reclass, template = z))

When calling both of these in the python script it is about an hour vs 40 min within Python, line by line in R is less than a minute a piece.

The main problem other than the time it takes is that my output raster, as mentioned, for the last 86 pixels contains NoData values.

which(is.na(Zstat))
results: integer(0)

HOWEVER, after the matrix as vector.

((reclass <- matrix(Zstat[,"vals"][as.vector(z)+1]))) 

is run. I now have NoData values where I should NOT……

head(which(is.na(reclass_2)))

>[1] 322595978 322596008 322596041 322596042 322596118 322596237

> length(reclass_2)
[1] 324036001
> na <- which(is.na(reclass_2))
> length(na)
[1] 32901

So what is causing 32901 pixels to be NA..?

Images:

First: right most green dot, in grey NoData area.

Second: left most green dot, in brown value zero area

third: dark dot in green area, green is NoData and should be value 1

Files mentioned:

Final_Zonal_Stat_file: R zonal output raster

ARCgis_Zonal output: created with ArcGIS model

Clump/Zones File: created in python with rpy2

water information: python gdal recoded raster file

My Zone NoData, Arc Zone output zero(white area, green dot)

My Zone NoData, Arc Zone output zero(brown area, green dot) My Zone NoData, Arc Zone output one (green area, dark dot)

UPDATE:

system.time(reclass <- matrix(Zstat[,"vals"][as.vector(z)+1]))

The clump file is missing sequential values for some reason, even though R was set to gaps=FALSE..... so my output for Zstat looks like this in the problem area:

         z    vals 
1872968 1872967 0
1872969 1872968 0
1872970 1872969 0
1872971 1884697 0 
1872972 1884698 0
1872973 1884699 0

The jump in numbers in row 1872971 is the problem. All the second column from this point on are not written to the final output raster. or reclass variable. and result in NoData values in my final matrix/raster.

Your Answer

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

Browse other questions tagged or ask your own question.