Skip to content

spgrass6, readVECT and overwriting ?

6 messages · Mathieu Rajerison, Rainer M Krug

Original
Rainer M Krug
# 
Mathieu Rajerison <mathieu.rajerison at gmail.com> writes:
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.
Don't worry abut your English - it is perfectly fine.
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.
OK - thanks.
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?
Sorry - overlooked it.
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

  
    

      
      
    

  


  


    

    
1 day later

    

      
      

        


  

        

      Mathieu Rajerison
    

    

      
      #
    

  


  

      
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
    

    

      
      #
    

  


  

      
Before was the code that works.

Below is the code that generates issues :

(sorry the multiple posts)

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"))

  cleaned = readVECT("clean")

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

2015-10-28 14:18 GMT+01:00 Mathieu Rajerison <mathieu.rajerison at gmail.com>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Rainer M Krug
    

    

      
      #
    

  


  

      
Mathieu Rajerison <mathieu.rajerison at gmail.com> writes:

            

      
      
      
    

        
OK. So what you are doing is in a loop

  write a raster to grass
  calculations in GRASS
  read result vector layer into R

And the error occurs in 

cleaned = readVECT("clean")

as temporary files during the readVECT are not cleaned up and not
overwritten - correct?

Then the rest is irrelevant.

You can reproduce it by doing

for (i in 1:10) {
  x <- readVECT("vectorLayerInSampleDataset")
}

Could you try this and see if it works in a sample dataset?

If it does, the error comes from somewhere else.


Under grass 7 it works:

,----
| 03:37:37 {master} ~/ownCompiled/emacs-mac$ grass70
| Cleaning up temporary files...
| Starting GRASS GIS...
| 
|           __________  ___   __________    _______________
|          / ____/ __ \/   | / ___/ ___/   / ____/  _/ ___/
|         / / __/ /_/ / /| | \__ \\_  \   / / __ / / \__ \
|        / /_/ / _, _/ ___ |___/ /__/ /  / /_/ // / ___/ /
|        \____/_/ |_/_/  |_/____/____/   \____/___//____/
| 
| Welcome to GRASS GIS 7.0.1
| GRASS GIS homepage:                      http://grass.osgeo.org
| This version running through:            Bash Shell (/usr/local/bin/bash)
| Help is available with the command:      g.manual -i
| See the licence terms with:              g.version -c
| If required, restart the GUI with:       g.gui wxpython
| When ready to quit enter:                exit
| 
| Launching <wxpython> GUI in the background, please wait...
| GRASS 7.0.1 (nc_basic_spm_grass7):~/ownCompiled/emacs-mac > R
| 
| R version 3.2.2 (2015-08-14) -- "Fire Safety"
| Copyright (C) 2015 The R Foundation for Statistical Computing
| Platform: x86_64-apple-darwin14.5.0 (64-bit)
| 
| R is free software and comes with ABSOLUTELY NO WARRANTY.
| You are welcome to redistribute it under certain conditions.
| Type 'license()' or 'licence()' for distribution details.
| 
|   Natural language support but running in an English locale
| 
| R is a collaborative project with many contributors.
| Type 'contributors()' for more information and
| 'citation()' on how to cite R or R packages in publications.
| 
| Type 'demo()' for some demos, 'help()' for on-line help, or
| 'help.start()' for an HTML browser interface to help.
| Type 'q()' to quit R.
| 
| > library(rgrass7)
| Loading required package: sp
| Loading required package: XML
| GRASS GIS interface loaded with GRASS version: GRASS 7.0.1 (2015)
| and location: nc_basic_spm_grass7
| > x <- readVECT("streams")
| Exporting 8554 features...
|  100%
| v.out.ogr complete. 8554 features (Line String type) written to <streams>
| (ESRI_Shapefile format).
| OGR data source with driver: ESRI Shapefile
| Source: "/Users/rainerkrug/nc_basic_spm_grass7/PERMANENT/.tmp/Rainers-MacBook-Pro.local", layer: "streams"
| with 8554 features
| It has 14 fields
| > x <- readVECT("streams")
| Exporting 8554 features...
|  100%
| v.out.ogr complete. 8554 features (Line String type) written to <streams>
| (ESRI_Shapefile format).
| OGR data source with driver: ESRI Shapefile
| Source: "/Users/rainerkrug/nc_basic_spm_grass7/PERMANENT/.tmp/Rainers-MacBook-Pro.local", layer: "streams"
| with 8554 features
| It has 14 fields
| >
`----


and also under grass-64:

,----
| > library(spgrass6)
| Loading required package: sp
| Loading required package: XML
| GRASS GIS interface loaded with GRASS version: GRASS 6.4.4 (2014)
| and location: spearfish60
| > x_readVECT("soils")
| Error: could not find function "x_readVECT"
| > x <- readVECT("soils")
| with_c: argument reversed from version 0.7-11 and in GRASS 6
| Exporting 737 areas (may take some time)...
|  100%
| WARNING: 1 features without attributes were written
| v.out.ogr complete. 737 features written to <soils> (ESRI_Shapefile).
| OGR data source with driver: ESRI Shapefile
| Source: "/Users/rainerkrug/spearfish60/user1/.tmp/Rainers-MacBook-Pro.local", layer: "soils"
| with 737 features
| It has 2 fields
| Warning message:
| In execGRASS("v.out.ogr", flags = flags, input = vname, type = type,  :
|   The command:
| v.out.ogr -e -c input=soils type=area layer=1 dsn=/Users/rainerkrug/spearfish60/user1/.tmp/Rainers-MacBook-Pro.local olayer=soils format=ESRI_Shapefile
| produced at least one warning during execution:
| Exporting 737 areas (may take some time)...
|  100%
| WARNING: 1 features without attributes were written
| v.out.ogr complete. 737 features written to <soils> (ESRI_Shapefile).
| > plot(x)
| > plot(x)
| > x <- readVECT("soils")
| with_c: argument reversed from version 0.7-11 and in GRASS 6
| Exporting 737 areas (may take some time)...
|  100%
| WARNING: 1 features without attributes were written
| v.out.ogr complete. 737 features written to <soils> (ESRI_Shapefile).
| OGR data source with driver: ESRI Shapefile
| Source: "/Users/rainerkrug/spearfish60/user1/.tmp/Rainers-MacBook-Pro.local", layer: "soils"
| with 737 features
| It has 2 fields
| Warning message:
| In execGRASS("v.out.ogr", flags = flags, input = vname, type = type,  :
|   The command:
| v.out.ogr -e -c input=soils type=area layer=1 dsn=/Users/rainerkrug/spearfish60/user1/.tmp/Rainers-MacBook-Pro.local olayer=soils format=ESRI_Shapefile
| produced at least one warning during execution:
| Exporting 737 areas (may take some time)...
|  100%
| WARNING: 1 features without attributes were written
| v.out.ogr complete. 737 features written to <soils> (ESRI_Shapefile).
| >
`----


Also: if you want to export it to a shape file, why do you import it
into R and save it afterwards? Why not write it directly from GRASS as a
shape file?

Cheers,

Rainer

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Mathieu Rajerison
    

    

      
      #
    

  


  

      
Yes, you're right. I should export directly from GRASS.

When doing :

for (i in 1:3) {
  f <- readVECT("f")
}

I get :
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERREUR :Layer <f> 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: "f"
with 1 features
It has 5 fields
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERREUR :Layer <f> 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: "f"
with 1 features
It has 5 fields
with_c: argument reversed from version 0.7-11 and in GRASS 6
ERREUR :Layer <f> 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: "f"
with 1 features
It has 5 fields
Warning messages:
1: running command 'v.out.ogr.exe -e -c input=f type=area layer=1
dsn=D:/GRASSDB/paca/mapset/.tmp olayer=f format=ESRI_Shapefile' had status
1
2: running command 'v.out.ogr.exe -e -c input=f type=area layer=1
dsn=D:/GRASSDB/paca/mapset/.tmp olayer=f format=ESRI_Shapefile' had status
1
3: running command 'v.out.ogr.exe -e -c input=f type=area layer=1
dsn=D:/GRASSDB/paca/mapset/.tmp olayer=f format=ESRI_Shapefile' had status
1

-----------------
sessionInfo() :
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252
 LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C                   LC_TIME=French_France.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] R.utils_2.1.0     R.oo_1.19.0       R.methodsS3_1.7.0 igraph_1.0.1
 gdalUtils_2.0.1.7
 [6] spdep_0.5-88      Matrix_1.2-2      maptools_0.8-34   spgrass6_0.8-8
 XML_3.98-1.3
[11] rgeos_0.3-8       FNN_1.1           rgdal_0.9-2       raster_2.3-40
  sp_1.0-17
[16] RStoolbox_0.1.1

loaded via a namespace (and not attached):
 [1] boot_1.3-13        car_2.0-25         caret_6.0-57       coda_0.18-1
     codetools_0.2-9
 [6] colorspace_1.2-6   deldir_0.1-9       digest_0.6.8
doParallel_1.0.8   foreach_1.4.2
[11] foreign_0.8-61     geosphere_1.4-3    ggplot2_1.0.1      grid_3.1.2
      gtable_0.1.2
[16] iterators_1.0.7    lattice_0.20-29    LearnBayes_2.15    lme4_1.1-10
     magrittr_1.5
[21] MASS_7.3-35        MatrixModels_0.4-1 mgcv_1.8-3         minqa_1.2.4
     munsell_0.4.2
[26] nlme_3.1-118       nloptr_1.0.4       nnet_7.3-8
parallel_3.1.2     pbkrtest_0.4-2
[31] plyr_1.8.3         proto_0.3-10       quantreg_5.19      Rcpp_0.12.0
     reshape2_1.4.1
[36] scales_0.2.5       SparseM_1.7        splines_3.1.2      stats4_3.1.2
      stringi_0.5-5
[41] stringr_1.0.0      tools_3.1.2

2015-10-28 16:16 GMT+01:00 Rainer M Krug <Rainer at krugs.de>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Rainer M Krug
    

    

      
      #
    

  


  

      
Mathieu Rajerison <mathieu.rajerison at gmail.com> writes:

            

      
      
      
    

        
This looks wrong. Could you please run a reproducible example using one
of the sample datasets, i.e.

1) start grass 6.4 and selsct a sample dataset as kocation
2) in the grass session start R
3) library(spgrass6)
4) for (i in 1:3) {x <- readVECT("soil")}

and send the transcript and the error messages?

Thanks,

Rainer