Skip to content
Prev 23608 / 29559 Next

spgrass6, readVECT and overwriting ?

Mr Krug,


It is difficult to reproduce an example from the spearfish dataset as it
doesn't natively contain data in it close to the one I used.

To give you th context, my goal is to detect all blue features : rivers
from a topographic map, ie 3-band RGB raster.

With RStoolbox, I perform a Spectral Angle Classification of my data using
sample points located on (near)blue pixels

With GRASS, I then thin, vectorize and clean the result.

As the process must be performed on a very large scale, I split it using
grids where sample RGB parts of the topo map are extracted with
gdal_translate. Everything is executed within a loop. At the end of the
loop, I export the resulted cleant shapefile in a folder. At the end, I
have as many shapefiles as processing grids.

That's why I use temporary vector datasets inside GRASS that are
overwritten at each loop instead of storing each individual result. I
noticed that in my process, the vector data read in the grass database was
not correct : always the first one. Because in fact, it didn't manage to
delete the temporary file before writing the new one..

The best I can do, is give you part of the code as as you can see the
context.

Thanks for the code improvements with file.remove and unlink !

Cheers,

Mathieu

for (id in 1:length(grids)) {

(...)

## SMOOTH
  print(">> SMOOTH")
  smoothed = smooth(sieved, 9)
  values(smoothed)[values(smoothed) > 0] = 1
  values(smoothed)[values(smoothed)==0] = NA

  ## THIN
  # WRITE TO GRASS
  writeRAST(as(smoothed, "SpatialGridDataFrame"), "bdtopo", zcol="layer",
flags=c("overwrite"))
  execGRASS("g.region", rast="bdtopo")
  execGRASS("r.thin", parameters=list(input="bdtopo", output="thin"),
flags=c("overwrite"))
  thin = raster(readRAST("thin"))

  ## CLEAN
  execGRASS("r.to.vect", parameters=list(input="thin", output="thin",
feature="line"), flags=c("overwrite"))
  execGRASS("v.clean", parameters=list(input="thin", tool="rmdangle",
output="clean", thresh=20), flags=c("overwrite"))

  # REMOVE FILES
  G = gmeta()
  tmpDir = file.path(G$GISDBASE, G$LOCATION_NAME, G$MAPSET, ".tmp")
  file.remove( list.files(path = tmpDir, Pattern ="^clean.*$", full.names =
TRUE))

  # READ TO SPDF
  execGRASS("v.out.ogr", parameters=list(input="clean", type="line",
dsn=tmpDir, olayer="clean", format="ESRI_Shapefile"))
  cleaned = readOGR(tmpDir, "clean")

  ## EXPORT
  writeOGR(cleaned, "OUT/CLEANED", id, "ESRI Shapefile")
}

2015-10-27 14:11 GMT+01:00 Rainer M Krug <Rainer at krugs.de>: