I've decided to process my Landsat data in R instead or ArcGIS - due to my missing knowledge of python and and because of (assumed) high computation capacities of R. I want to :
- import r1 raster to R,
- import shp1 convert raster r1 to shp
r.to.poly (dissolve = TRUE)
- intersect converter raster
r.to.polywith my polygon shp1
- calculate area of every created polygon of intersected shp
# read shp shp <-readOGR(dsn = "C://...", layer = "m") #read raster r1<-raster("r1.tif") # convert raster to polygon, dissolved neighboring same values r.to.poly<-rasterToPolygons(r1, dissolve = T) # define the same projection proj4string(shp) <- proj4string(r.to.poly) # use intersection from raster package int.r <-raster::intersect(r.to.poly,shp) # calculate area per polygon int.r$area <-gArea(int.r, byid = T) # export shapefile writeOGR(int.r, dsn = "C:/...", layer = "...", driver="ESRI Shapefile", overwrite = TRUE)
So far, so good, but I takes about an hour to run the single conversion ! moreover, when I tried FOR loop, my R on Windows crashed twice... It runs on mac, for the moment. Where the problem could be and how can I increase computation speed? Am I running out of R memory? The raster size on my disk is only 779 580 byte, size of shp is 1 729 532 bytes, thus are small. Also, make the same task in ArcGIS takes only couple seconds.
I've found some related discussion here: Increasing speed of crop, mask, & extract raster by many polygons in R? but as I have only about 10 rasters to process I don't want to start with parallel processing...