Saturday 23 February 2019

python - Efficiently reading and reclassifying many rasters in R?


I've been tasked to create a suitability analysis of wave conditions in the Gulf of Mexico. I have 2 thousand or so raster files that are about 8 MB each (2438 columns, 1749 rows, 1km cell size). The parameter in the raster files is wave period and I'd like to reclassify all of the rasters such that if 4<= wave period <=9 then make cell = 1, else cell = 0. Then sum up all the rasters into a final raster and divide by the total number of rasters to give a total percentage of suitable observations and export result into some ESRI-compatible format...maybe something that can support floats if need be. I've not worked much at all with either Python or R, but after searching online it seems to make sense to perform this process in one of those languages. I've come up with some code so far in R, but am confused about how to make this work.



library(rgdal)
raster_data <- list.files(path=getwd()) #promt user for dir containing raster files
num_files <- length(raster_data)
for (i in raster_data) { #read in rasters
my_data <- readGDAL(raster_data[i])

At this point I'm confused as to whether I should also reclassify and start summing the data within this loop or not. My guess would be yes since otherwise I think I would possibly run out of memory on my computer, but just not sure. I'm also not sure about how to reclassify the data.


In researching online I think I would use reclass(my_data, c(-Inf,3,0, 4,9,1, 10,Inf,0)), but does that look right?


And for summing would I use sum(stack(my_data)) and somehow output that. Also...if this might be more efficiently performed or written in Python I'd be open to that as well. I truly am a beginner when it comes to programming.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...