[R-sig-Geo] spgrass6, readVECT and overwriting ?

Mathieu Rajerison mathieu.rajerison at gmail.com
Wed Oct 28 14:18:20 CET 2015 

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>:

> Mathieu Rajerison <mathieu.rajerison at gmail.com> writes:
>
> > Dear Mr Krug,
> >
>
> I'll Cc the R-sig-geo list in to keep the conversation there for future
> reference and for other users who might have the same problem.
>
> > Sorry in advance for my english.
>
> Don't worry abut your English - it is perfectly fine.
>
> >
> > My question was not about debugging, but about a potentiel generic
> > behaviour of spgrass6.
> >
> > I should have reformulated my question in a simpler manner : does
> spgrass6
> > allow overwriting an existing dataset when using readVECT function ?
>
> The functions readVect / RAST are using temporary files - you are right,
> but they should delete these upon exit. So if they do nbot do this, it
> is a bug.
>
> >
> > 1) no version information about R and GRASS (there are several GRASS 6
> > versions)
> > I've just put it at the end of this post.
> > 2) no session information i.e. package versions used
> > the latest spgrass6
> > 3) no information about the OS (I guess it is windows based on the
> > Source path above)
> > you're right
>
> OK - thanks.
>
> > 4) no reproducible example which could be easily done by using the GRASS
>
> As I said - if you manually have to delete the temporary files, then
> there is a bug in the package. Could you please provide a reproducible
> example, using one of the example data sets from GRASS, so that I can
> look at it?
>
> > 5) as far as I can see, in the newest version of spgrass6 there is no
> > function readVECT()
> > Yes, there is, as you will see on page 14 of the documentation.
>
> Sorry - overlooked it.
>
> >
> > Finally, I solved the problem by using the following hint, "clean" being
> > the name of my temporary shapefile dataset :
> >   l = list.files("D:/GRASSDB/paca/mapset/.tmp", "^clean.*$"); l =
> > file.path("D:/GRASSDB/paca/mapset/.tmp", l)
> >   lapply(l,  file.remove)
> >   execGRASS("v.out.ogr", parameters=list(input="clean", type="line",
> > dsn="D:/GRASSDB/paca/mapset/.tmp", olayer="clean",
> format="ESRI_Shapefile"))
>
> OK - reformated, the code looks as follow:
>
> --8<---------------cut here---------------start------------->8---
> l = list.files("D:/GRASSDB/paca/mapset/.tmp", "^clean.*$")
> l = file.path("D:/GRASSDB/paca/mapset/.tmp", l)
> lapply(l,  file.remove)
> execGRASS( "v.out.ogr",
>            parameters=list(input="clean",
> type="line",dsn="D:/GRASSDB/paca/mapset/.tmp", olayer="clean",
> format="ESRI_Shapefile"))
> --8<---------------cut here---------------end--------------->8---
>
> as a side note, you should be able to do
>
> file.remove( list.files(path = "D:/GRASSDB/paca/mapset/.tmp", Pattern =
> "^clean.*$", full.names = TRUE) )
>
> or
>
> unlink("D:/GRASSDB/paca/mapset/.tmp/clean.*")
>
> instead of the first three commands.
>
> Now - in your example, you don't use the readVECT() function?
>
> If you check in the GRASS help for v.out.ogr, you can specify the
> overwrite flag and the exported layer will be overwritten.
>
> But I have the feeling you are doing something in a to complicated way.
>
> What is it you want to do, and e=what is the problem you have or the
> error message you get?
>
> Cheers,
>
> Rainer
>
> >
> >
> >> R.Version()
> > $platform
> > [1] "x86_64-w64-mingw32"
> > $arch
> > [1] "x86_64"
> > $os
> > [1] "mingw32"
> >
> > $system
> > [1] "x86_64, mingw32"
> > (...)
> > $version.string
> > [1] "R version 3.1.2 (2014-10-31)"
> > $nickname
> > [1] "Pumpkin Helmet"
> >
> > I will happily look int this if you can provide the necessary info.
> >
> > Cheers,
> >
> > Mathieu
> >
> > 2015-10-24 13:34 GMT+02:00 Rainer M Krug <Rainer at krugs.de>:
> >
> >> Mathieu Rajerison <mathieu.rajerison at gmail.com> writes:
> >>
> >> > Hi,
> >> >
> >> > I'm using spgrass6 and I use readVECT function in a loop.
> >> >
> >> > Using it causes an error because it doesn't overwrite the the
> temporary
> >> > shapefile..
> >> >
> >> >> thin = readVECT(vname="zoneBDTOPOthin")
> >> > with_c: argument reversed from version 0.7-11 and in GRASS 6
> >> > ERROR :Layer <zoneBDTO> already exists in OGR data source
> >> >         'D:/GRASSDB/paca/mapset/.tmp'
> >> > OGR data source with driver: ESRI Shapefile
> >> > Source: "D:/GRASSDB/paca/mapset/.tmp", layer: "zoneBDTO"
> >> > with 770 features
> >> > It has 3 fields
> >> > Warning message:
> >> > running command 'v.out.ogr.exe -e -c input=zoneBDTOPOthin2 type=line
> >> > layer=1 dsn=D:/GRASSDB/paca/mapset/.tmp olayer=zoneBDTO
> >> > format=ESRI_Shapefile' had status 1
> >> >
> >> >
> >> > I wonder how to specify to spgrass that I want to overwrite the file
> in
> >> > readVECT, although I know using execGRASS and v.out.ogr with the
> >> > appropriate overwrite flags would do the job.
> >> >
> >> > Thanks in advance for the answer..
> >>
> >> One reason why probably nobody replied is that
> >>
> >> 1) no version information about R and GRASS (there are several GRASS 6
> >> versions)
> >> 2) no session information i.e. package versions used
> >> 3) no information about the OS (I guess it is windows based on the
> >> Source path above)
> >> 4) no reproducible example which could be easily done by using the GRASS
> >> 5) as far as I can see, in the newest version of spgrass6 there is no
> >> function readVECT()
> >>
> >> So I would suggest to provide the info and give a reproducible example
> >> using a GRASS sample dataset
> >> (see https://grass.osgeo.org/download/sample-data/ - the Spearfish
> >> dataset is probably the best)
> >>
> >> I will happily look int this if you can provide the necessary info.
> >>
> >> Cheers,
> >>
> >> Rainer
> >>
> >> >
> >> > Mathieu
> >> >
> >> >       [[alternative HTML version deleted]]
> >> >
> >> > _______________________________________________
> >> > R-sig-Geo mailing list
> >> > R-sig-Geo at r-project.org
> >> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >> --
> >> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
> >> Biology, UCT), Dipl. Phys. (Germany)
> >>
> >> Centre of Excellence for Invasion Biology
> >> Stellenbosch University
> >> South Africa
> >>
> >> Tel :       +33 - (0)9 53 10 27 44
> >> Cell:       +33 - (0)6 85 62 59 98
> >> Fax :       +33 - (0)9 58 10 27 44
> >>
> >> Fax (D):    +49 - (0)3 21 21 25 22 44
> >>
> >> email:      Rainer at krugs.de
> >>
> >> Skype:      RMkrug
> >>
> >> PGP: 0x0F52F982
> >>
>
> --
> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
> Biology, UCT), Dipl. Phys. (Germany)
>
> Centre of Excellence for Invasion Biology
> Stellenbosch University
> South Africa
>
> Tel :       +33 - (0)9 53 10 27 44
> Cell:       +33 - (0)6 85 62 59 98
> Fax :       +33 - (0)9 58 10 27 44
>
> Fax (D):    +49 - (0)3 21 21 25 22 44
>
> email:      Rainer at krugs.de
>
> Skype:      RMkrug
>
> PGP: 0x0F52F982
>

	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list