7

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 :

  1. import r1 raster to R,
  2. import shp1 convert raster r1 to shp r.to.poly (dissolve = TRUE)
  3. intersect converter raster r.to.poly with my polygon shp1
  4. calculate area of every created polygon of intersected shp

Thus:

# 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...

  • If the speed isn't acceptable perhaps try GDAL_Polygonize gdal.org/gdal_polygonize.html as a process.. or as you have access to ArcGis but lack experience in arcpy perhaps I could help with that instead. – Michael Stimson May 8 '16 at 23:31
  • thank you @MichaelMiles-Stimson, I've solved in with ModelBuilder for a while, it was just 11 rasters. But I will really appreciate your help in the future ! ;) – maycca May 9 '16 at 18:06
  • 1
    Which step is taking the most time? Converting to polygons, or the intersection? – amball May 23 '16 at 4:47
  • 1
    Also, are you only interested in the part of your raster that intersects with your shapefile? If so, you might try using raster::intersect to crop your raster BEFORE converting it to polygons. – amball May 23 '16 at 4:55
5
+50

Gdal Installation

Install Gdal command line tools and check to see if its binaries are added to path environment variable. e.g. in windows: open Run and type:

rundll32.exe sysdm.cpl,EditEnvironmentVariables 

Then follow the screenshot

enter image description here

Download and install gdal python bindings from here according to your python and OS.

install it using:

pip.exe install GDAL-2.0.2-cp27-none-win32.whl

You may encounter issues while installing gdal. Please see Installing GDAL with Python on windows?

In that thread users have suggested that Gdal binaries from gisinternals installs both the command line tools and the python bindings. Try to install it from there. Thus none of the above steps would be relevant.

To Make sure Gdal is installed open a command prompt and type:

ogrinfo

and check it works and you don't get:

'ogrinfo' is not recognized as an internal or external command, operable program or batch file.

To check gdal python bindings is installed open command prompt and type:

python
import gdal

if you get the following error then it is not installed:

Traceback (most recent call last): File "", line 1, in ImportError: No module named gdal

R side

Source the following r function from John Baumgartner blog:

gdal_polygonizeR <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile',
                             pypath=NULL, readpoly=TRUE, quiet=TRUE) {
  if (isTRUE(readpoly)) require(rgdal)
  if (is.null(pypath)) {
    pypath <- Sys.which('gdal_polygonize.py')
  }
  if (!file.exists(pypath)) stop("Can't find gdal_polygonize.py on your system.")
  owd <- getwd()
  on.exit(setwd(owd))
  setwd(dirname(pypath))
  if (!is.null(outshape)) {
    outshape <- sub('\\.shp$', '', outshape)
    f.exists <- file.exists(paste(outshape, c('shp', 'shx', 'dbf'), sep='.'))
    if (any(f.exists))
      stop(sprintf('File already exists: %s',
                   toString(paste(outshape, c('shp', 'shx', 'dbf'),
                                  sep='.')[f.exists])), call.=FALSE)
  } else outshape <- tempfile()
  if (is(x, 'Raster')) {
    require(raster)
    writeRaster(x, {f <- tempfile(fileext='.tif')})
    rastpath <- normalizePath(f)
  } else if (is.character(x)) {
    rastpath <- normalizePath(x)
  } else stop('x must be a file path (character string), or a Raster object.')
  system2('python', args=(sprintf('"%1$s" "%2$s" -f "%3$s" "%4$s.shp"',
                                  pypath, rastpath, gdalformat, outshape)))
  if (isTRUE(readpoly)) {
    shp <- readOGR(dirname(outshape), layer = basename(outshape), verbose=!quiet)
    return(shp)
  }
  return(NULL)
}

use this command to do the job:

r.to.poly<-gdal_polygonizeR(r1,pypath = "C:\\Program Files\\GDAL\\gdal_polygonize.py")#, dissolve = T)

Change pypath parameter according to your system. Be cautious that gdal_plygonize creates huge shapefiles. My 1 MB tif converted to a 128 MB shapefile. R needs a lot of memory to open this shapefile. Although the conversion was very fast. Thanks to python and gdal!

Another option would be Esri r-bridge to do the computation in Arcgis and return the output to R. However r-bridge doesn't support raster layers yet. (Thanks to @JeffreyEvans)

However the Gdal method is commercial free!

  • Have you used the ESRI R-bridge? In previous versions it was a one-way street from ArcGIS and not designed to call ArcGIS processes from R. – Jeffrey Evans May 23 '16 at 18:51
  • no, I didn't use R-bridge. canI use it just like an another library? and, moreover, I can avoid all of the GDAL installation? moreover, I am sure that i have python on my PC, as I hve GRASS and QGIS... – maycca May 23 '16 at 20:20
  • I'm sorry @FaridCher.. I can't follow it! where can I put 'path = %path%;C:\Program Files\GDAL'? in R? I/ve downloaded the GDAL-2.0.2-cp35-none-win32.whl, and also I've installed OSGeo4W, but I am sure that 've already had this because I am using QGIS... also, when I tried to install pip.exe, I couldn't find the file ! please, is your code chronological, or there is different order of steps? thank you ! – maycca May 23 '16 at 20:30
  • @JeffreyEvans I haven't tried it myself just a suggestion on top of my head. I checked the latest version yet native arcgis raster datasets are not supported. However one could use Arcpy and create a python script like gdal_polygonizeR.py and then call it from r and use the supported raster formats by R as the output. – Farid Cheraghi May 23 '16 at 22:58
  • 1
    Thank you a lot @FaridCher, it works well now!! My problem was, that I was ustaling it while RStudio was running, and I was testing it on R. After Closing all programs and assigning path to MAC bash: echo 'export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH' >> ~/.bash_profile and defining my pypath = "/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py" after sourcing the gdal_polygonizeR it worked well. Thank you again for all help ! My 30 minutes needed to convert raster to polygon takes now only 3 seconds !! ;) – maycca May 27 '16 at 18:19
0

I used rasterToPolygons from the raster package too in the past, but now I prefer gdal_polygonizeR by John Baumgartner. It bases on gdal_polygonize.py and is much faster. John Baumgartner published the code and gave an example for usage in his blog.

If you are familiar with python you could use gdal_polygonize.py directly of course.

0

from : https://sites.google.com/site/eospansite/installinggdalonmacosx


I'm copying these things in case link be broken:

Installing GDAL on Mac OSX and Polygonize Rasters

GDAL Complete: you must install both gdal.pkg and numpy.pkg

In all cases:

1.1. Double click and you will get a *.dmg file in Downloads (or Descargas). Wait until the download completes (GDAL Complete will take longer).

1.2. Once the dmg is downloaded, it normally gets uncompressed by itself and the folder is automatically opened. Otherwise, go to Downloads and double-click on the dmg file.

1.3. Look for the *.pkg inside the folder that opens, double click, accept the license and proceed. Sometimes you have to click on the path for the installation in order to get the Accept button active.

Once you have all 4 packages installed,

  1. Open a Terminal (Applications/Utilities/Terminal) and write export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH and hit Return


READY TO GO TO R! :)

Read libraries:

library(raster)
library(rgeos)

Define the function aavailable from https://johnbaumgartner.wordpress.com/2012/07/26/getting-rasters-into-shape-from-r/

gdal_polygonizeR <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile',
                             pypath=NULL, readpoly=TRUE, quiet=TRUE) {
  if (isTRUE(readpoly)) require(rgdal)
  if (is.null(pypath)) {
    pypath <- Sys.which('gdal_polygonize.py')
  }
  if (!file.exists(pypath)) stop("Can't find gdal_polygonize.py on your system.")
  owd <- getwd()
  on.exit(setwd(owd))
  setwd(dirname(pypath))
  if (!is.null(outshape)) {
    outshape <- sub('\\.shp$', '', outshape)
    f.exists <- file.exists(paste(outshape, c('shp', 'shx', 'dbf'), sep='.'))
    if (any(f.exists))
      stop(sprintf('File already exists: %s',
                   toString(paste(outshape, c('shp', 'shx', 'dbf'),
                                  sep='.')[f.exists])), call.=FALSE)
  } else outshape <- tempfile()
  if (is(x, 'Raster')) {
    require(raster)
    writeRaster(x, {f <- tempfile(fileext='.tif')})
    rastpath <- normalizePath(f)
  } else if (is.character(x)) {
    rastpath <- normalizePath(x)
  } else stop('x must be a file path (character string), or a Raster object.')
  system2('python', args=(sprintf('"%1$s" "%2$s" -f "%3$s" "%4$s.shp"',
                                  pypath, rastpath, gdalformat, outshape)))
  if (isTRUE(readpoly)) {
    shp <- readOGR(dirname(outshape), layer = basename(outshape), verbose=!quiet)
    return(shp)
  }
  return(NULL)
}

make raster polygonization on dummy raster

- the way you define the path may be different on R, RStudio and R from Terminal !! (at least on OSX 10.9.5)

# read raster
r5<-raster("r5.img")
plot(r5)


# in R + RStudio: define whole path where to find gdal_polygonize.py !!
r.to.poly<-gdal_polygonizeR(r5,pypath = "/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py")#


#in R run from Terminal:  no need to define the whole path !!
r.to.poly<-gdal_polygonizeR(r5)

to check, if you need to define the path to gdal_polygonize.py or not:

> Sys.which("gdal_polygonize.py")
                                              gdal_polygonize.py 
"/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py" # means path is defined, you should be able just to run 

r.to.poly<-gdal_polygonizeR(r5)

If it looks like this:

> Sys.which("gdal_polygonize.py")
gdal_polygonize.py 
                "" 

you would need to define the whole path to access gdal_polygonize.py

r.to.poly<-gdal_polygonizeR(r5,pypath = "/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py")
0

You can do the analysis described in your post without converting the raster to a polygon. Use the raster::extract function to extract the raster values to each polygon. You can then use lapply on the resulting list object with table to return cell counts of each class. For area of each raster class, you just use a standard conversion of cell area and counts. Here is a quick example.

First lets create a discrete raster with values [5:8] and two polygons.

library(raster)

r <- raster(ncol=36, nrow=18)
  r[] <- round(runif(ncell(r),5,8),0)

polys <- spPolygons(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)), 
                    rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)))

plot(r)
  plot(polys, add=TRUE)

Now we can extract the raster values for each polygon. The lapply function is used with table to get the cell counts of each raster class.

( v <- extract(r, polys) )
( v.counts <- lapply(v, FUN = table) )

We can now merge the raster cell counts with the polygon data. The do.call function is used to create a data.frame from the list object. The list is ordered the same as the polygons so it can be directly merged. The columns are cell counts for each raster class value intersecting a given polygon.

polys <- SpatialPolygonsDataFrame(polys, data.frame(IDS=1:length(polys), 
                                  do.call("rbind", v.counts)))
  names(polys@data)[2:ncol(polys)] <- paste("class", 5:8, sep=".")
  polys@data

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.