Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rast ng probably ready to merge #45

Merged
merged 20 commits into from
Mar 4, 2022
Merged

Rast ng probably ready to merge #45

merged 20 commits into from
Mar 4, 2022

Conversation

rsbivand
Copy link
Owner

@rsbivand rsbivand commented Feb 18, 2022

The rastNG branch provides re-implementations using the "SpatRaster" and "SpatVector" classes from terra, and using terra methods for reading and writing intermediate files (See #42). The additions need checking on (R CMD check rgrass7_0.2-8.tar.gz in nc_basic_spm_grass7 location):
rgrass7_0.2-8.tar.gz
or this example file:
rgrass7-Ex.zip

  • Windows standalone GRASS 8.0 binary R 4.1.2
  • Windows OSGeo4W GRASS 8.0 binary R 4.1.2
  • macOS x86_64 GRASS 8.0 binary (via Rosetta on M1) R 4.1.2
  • macOS arm64 (arm64 GRASS not available)
  • Windows standalone GRASS 7.8 binary R 4.1.2
  • Windows OSGeo4W GRASS 7.8 binary R 4.1.2
  • macOS x86_64 GRASS 7.8 binary (via Rosetta on M1) R 4.1.2
  • Fedora GRASS 8.0 R 4.1.2
  • Fedora GRASS 7.8 R 4.1.2
  • Fedora GRASS 8.0 R 4.2dev
  • Fedora GRASS 7.8 R 4.2dev

before merging (possibly other Linux versions if contributions are possible).

Copy link
Collaborator

@hellik hellik left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no chance to test it at the moment. though, it looks very promising.

Copy link
Collaborator

@veroandreo veroandreo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @rsbivand! This is great! I downloaded the rgrass7_0.2-8.tar.gz and ran R CMD check rgrass7_0.2-8.tar.gz in nc_basic_spm_grass7 location within a GRASS 8-dev session. All went smoothly :)

Just to clarify usage then (and check if I got it right):

The legacy use implying the use of use_sp() and use_sf() together with readRAST()/writeRAST() and readVECT()/writeVECT() (use_sf only for vectors though) remains the same for now but will be deprecated in the next version.

The new functionality in this PR adds terra support for both raster and vector formats with the functions read_VECT()/write_VECT() and read_RAST()/write_RAST(), and for raster, there's the option to set return_format="terra" or return_format="SGDF" (to keep sp support), but for vectors no specification is needed and by default terra's SpatVector class will be used. Did I understand it ok? I assume vector objects can then be coerced to sf if desired as you did for rasters and stars in the example in #42, right?

@rsbivand
Copy link
Owner Author

Yes, the path I'm suggesting is to deprecate read/write/RAST/VECT without underscore soon, and only use terra leaving r.in/out.bin for sp raster. Coercion to other representation formats would then be up to the user. That will be easier to maintain, I think, but welcome discussion. If for example stars is a better fit for temporal multiband rasters, that could be added later.

@rsbivand
Copy link
Owner Author

@VLucet and @florisvdh Since no negative review, I'll merge tomorrow to rgrass7 and rgrass, and rgrass7 to main shortly thereafter, unless I hear why this should not be done.

@florisvdh
Copy link
Collaborator

@rsbivand sorry, I had missed your review request. I can look at it later today.

Copy link
Collaborator

@florisvdh florisvdh left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In current state (aa1cc02), R CMD check runs OK in GRASS 8.0.0RC2 and 7.8.7RC1 (Linux Mint 20) only if GRASS_INSTALLATION is set (due to today's aa1cc02), otherwise errors in both cases:

* checking examples ... ERROR
Running examples inrgrass7-Ex.Rfailed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: initGRASS
> ### Title: Initiate GRASS session
> ### Aliases: initGRASS get.GIS_LOCK set.GIS_LOCK unset.GIS_LOCK
> ###   unlink_.gislock remove_GISRC
> ### Keywords: spatial
> 
> ### ** Examples
> 
> GRASS_INSTALLATION <- Sys.getenv("GRASS_INSTALLATION")
> run <- FALSE
> if (nzchar(GRASS_INSTALLATION)) run <- file.info(GRASS_INSTALLATION)$isdir
> run <- require(terra, quietly=TRUE)
terra 1.5.12
> if (run) {
+ f <- system.file("ex/elev.tif", package="terra")
+ r <- rast(f)
+ plot(r, col=grDevices::terrain.colors(50))
+ }
> if (run) {
+ (loc <- initGRASS(GRASS_INSTALLATION, home=tempdir(), SG=r, override=TRUE))
+ }
Error in initGRASS(GRASS_INSTALLATION, home = tempdir(), SG = r, override = TRUE) : 
   not found
Execution halted

It's nice that terra can serve as the common tool for raster and vector usage in rgrass/rgrass7. rgrass7 in this branch works with both GRASS 7 and GRASS 8, cf. below.

Output from code in #42 using GRASS 7.8.7RC1
> # after https://github.com/rsbivand/rgrass/issues/42#issuecomment-1040805090
> 
> library(rgrass7)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: (GRASS not running)
> my_GRASS <- "/usr/local/grass78"
> loc <- initGRASS(my_GRASS, home=tempdir(), gisDbase=".", 
+                  location="nc_basic_spm_grass7", mapset="user1", override=TRUE)
> 
> execGRASS("g.version")
GRASS 7.8.7RC1 (2022)
> execGRASS("g.gisenv")
GISDBASE='.';
LOCATION_NAME='nc_basic_spm_grass7';
MAPSET='user1';
GRASS_GUI='text';
> 
> execGRASS("g.region", raster="elevation")
> execGRASS("g.region", flags="p")
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228500
south:      215000
west:       630000
east:       645000
nsres:      10
ewres:      10
rows:       1350
cols:       1500
cells:      2025000
> 
> execGRASS("r.slope.aspect", flags="overwrite", elevation="elevation",
+           slope="slope", aspect="aspect")
 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete
> use_sp()
> u00 <- readRAST(c("elevation", "slope", "aspect")) # legacy
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Warning message:
In showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum NAD83 (High Accuracy Reference Network) in Proj4 definition
> object.size(u00) # in memory
48607248 bytes
> u0 <- rgrass7::read_RAST(c("elevation", "slope", "aspect"),
+                           return_format="SGDF") # re-written SGDF
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
> object.size(u0) # in memory
48607248 bytes
> library(sp)
> all.equal(u00[["elevation"]], u0[["elevation"]])
[1] TRUE
> u1 <- rgrass7:::read_RAST(c("elevation", "slope", "aspect"),
+                           return_format="terra") # re-written terra
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpmiKrAM/file846235417544.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpmiKrAM/file846276a689b2.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpmiKrAM/file84626bb1b22d.grd> created.
> object.size(u1) # on file
1304 bytes
> library(terra)
terra 1.5.12
> u1
class       : SpatRaster 
dimensions  : 1350, 1500, 3  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_0=33.75 +lon_0=-79 +lat_1=36.1666666666667 +lat_2=34.3333333333333 +x_0=609601.22 +y_0=0 +ellps=GRS80 +units=m +no_defs 
sources     : file846235417544.grd  
              file846276a689b2.grd  
              file84626bb1b22d.grd  
names       : elevation,     slope,    aspect 
min values  :  55.57879,   0.00000,   0.00000 
max values  : 156.32986,  38.68939, 360.00000 
> all.equal(u00[["elevation"]], c(values(u1[["elevation"]])))
[1] TRUE
> library(stars)
Loading required package: abind
Loading required package: sf
Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is TRUE
> u2 <- st_as_stars(u1) # coerce to stars
> u2
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                          Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
file846235417544.grd  55.57879 97.80552 118.1428 112.9484 130.4342 153.3789
dimension(s):
           from   to offset delta                       refsys point
x             1 1500 630000    10 +proj=lcc +lat_0=33.75 +l...    NA
y             1 1350 228500   -10 +proj=lcc +lat_0=33.75 +l...    NA
attributes    1    3     NA    NA                           NA    NA
                                    values x/y
x                                     NULL [x]
y                                     NULL [y]
attributes elevation, slope    , aspect       
> object.size(u2) # in memory
48611744 bytes
> all.equal(u00[["elevation"]], c(dplyr::pull(u2[,,,"elevation"])))
[1] TRUE
> u2a <- read_stars(u1@ptr$filenames, proxy=TRUE)
> u2a
stars_proxy object with 3 attributes in 3 file(s):
$file846235417544.grd
[1] "[...]/file846235417544.grd"

$file846276a689b2.grd
[1] "[...]/file846276a689b2.grd"

$file84626bb1b22d.grd
[1] "[...]/file84626bb1b22d.grd"

dimension(s):
  from   to offset delta                       refsys point values x/y
x    1 1500 630000    10 +proj=lcc +lat_0=33.75 +l...    NA   NULL [x]
y    1 1350 228500   -10 +proj=lcc +lat_0=33.75 +l...    NA   NULL [y]
> object.size(u2a) # on file
10328 bytes
> all.equal(u00[["elevation"]], c(st_as_stars(u2a)[[1]])) # moved to memory
[1] TRUE
> 
> 
> 
> # after https://github.com/rsbivand/rgrass/issues/42#issuecomment-1043311481
> 
> f <- system.file("ex/elev.tif", package="terra")
> r <- rast(f)
> plot(r, col=grDevices::terrain.colors(50))
> loc <- initGRASS(my_GRASS, home=tempdir(), SG=r, override=TRUE)
> 
> execGRASS("g.gisenv")
GISDBASE='/tmp/RtmpmiKrAM';
LOCATION_NAME='file846267e0fe1';
MAPSET='file8462649e6aac';
GRASS_GUI='text';
> 
> rgrass7:::write_RAST(r, "elev", flags="overwrite")
Importing raster map <elev>...
 100%
> execGRASS("r.info", map="elev")
 +----------------------------------------------------------------------------+
 | Map:      elev                           Date: Mon Feb 28 18:35:54 2022    |
 | Mapset:   file8462649e6aac               Login of Creator: floris          |
 | Location: file846267e0fe1                                                  |
 | DataBase: /tmp/RtmpmiKrAM                                                  |
 | Title:                                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    FCELL                                                      |
 |   Rows:         90                                                         |
 |   Columns:      95                                                         |
 |   Total Cells:  8550                                                       |
 |        Projection: Latitude-Longitude                                      |
 |            N:  50:11:30N    S:  49:26:30N   Res: 0:00:30                   |
 |            E:      6:32E    W:   5:44:30E   Res: 0:00:30                   |
 |   Range of data:    min = 141  max = 547                                   |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal --overwrite input="/tmp/RtmpmiKrAM/file8462679681ca.grd" o\   |
 |    utput="elev" memory=300 offset=0 num_digits=0                           |
 |                                                                            |
 +----------------------------------------------------------------------------+

> execGRASS("r.slope.aspect", flags="overwrite", elevation="elev", slope="slope", aspect="aspect")
 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete
> u1 <- rgrass7:::read_RAST(c("elev", "slope", "aspect"), return_format="terra")
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpmiKrAM/file84625311d345.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpmiKrAM/file84626c197306.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpmiKrAM/file8462488e4019.grd> created.
> plot(u1[["elev"]], col=grDevices::terrain.colors(50))
> 
> 
> 
> # after https://github.com/rsbivand/rgrass/issues/42#issuecomment-1044769488
> 
> v1 <- terra::vect(system.file("ex/lux.shp", package = "terra"))
> v1
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 12, 6  (geometries, attributes)
 extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :  ID_1   NAME_1  ID_2   NAME_2  AREA       POP
 type        : <num>    <chr> <num>    <chr> <num>     <num>
 values      :     1 Diekirch     1 Clervaux   312 1.808e+04
                   1 Diekirch     2 Diekirch   218 3.254e+04
                   1 Diekirch     3  Redange   259 1.866e+04
> rgrass7:::write_VECT(v1, "lux")
Check if OGR layer <file84626201ccfa> contains polygons...
 100%
Creating attribute table for layer <file84626201ccfa>...
Default driver / database set to:
driver: sqlite
database: $GISDBASE/$LOCATION_NAME/$MAPSET/sqlite/sqlite.db
Importing 12 features (OGR layer <file84626201ccfa>)...
 100%
-----------------------------------------------------
Registering primitives...

-----------------------------------------------------
Cleaning polygons
-----------------------------------------------------
Breaking polygons...
Breaking polygons (pass 1: select break points)...
 100%
Breaking polygons (pass 2: break at selected points)...
 100%
-----------------------------------------------------
Removing duplicates...
 100%
-----------------------------------------------------
Breaking boundaries...
 100%
-----------------------------------------------------
Removing duplicates...
 100%
-----------------------------------------------------
Cleaning boundaries at nodes...
 100%
-----------------------------------------------------
Merging boundaries...
 100%
-----------------------------------------------------
Changing boundary dangles to lines...
 100%
-----------------------------------------------------
Building areas...
-----------------------------------------------------
Changing boundary bridges to lines...
 100%
-----------------------------------------------------
Registering primitives...

Building areas...
 100%
Attaching islands...
 100%
-----------------------------------------------------
Finding centroids for OGR layer <file84626201ccfa>...
 100%
-----------------------------------------------------
Writing centroids...
 100%
-----------------------------------------------------
12 input polygons
Total area: 2.56481E+09 (12 areas)
-----------------------------------------------------
Copying features...
 100%
Building topology for vector map <lux@file8462649e6aac>...
Registering primitives...

Building areas...
 100%
Attaching islands...
 100%
Attaching centroids...
 100%
Warning message:
In x@ptr$write(filename, layer, filetype, overwrite[1], options) :
  GDAL Message 6: dataset /tmp/RtmpmiKrAM/file84626201ccfa.gpkg does not support layer creation option ENCODING
> execGRASS("v.info", map="lux")
 +----------------------------------------------------------------------------+
 | Name:            lux                                                       |
 | Mapset:          file8462649e6aac                                          |
 | Location:        file846267e0fe1                                           |
 | Database:        /tmp/RtmpmiKrAM                                           |
 | Title:                                                                     |
 | Map scale:       1:1                                                       |
 | Name of creator: floris                                                    |
 | Organization:                                                              |
 | Source date:     Mon Feb 28 18:35:55 2022                                  |
 | Timestamp (first layer): none                                              |
 |----------------------------------------------------------------------------|
 | Map format:      native                                                    |
 |----------------------------------------------------------------------------|
 |   Type of map: vector (level: 2)                                           |
 |                                                                            |
 |   Number of points:       0               Number of centroids:  12         |
 |   Number of lines:        0               Number of boundaries: 33         |
 |   Number of areas:        12              Number of islands:    1          |
 |                                                                            |
 |   Map is 3D:              No                                               |
 |   Number of dblinks:      1                                                |
 |                                                                            |
 |   Projection: Latitude-Longitude                                           |
 |                                                                            |
 |               N:   50:10:53.83758N    S:  49:26:52.106316N                 |
 |               E:   6:31:41.707632E    W:    5:44:38.90454E                 |
 |                                                                            |
 |   Digitization threshold: 0                                                |
 |   Comment:                                                                 |
 |                                                                            |
 +----------------------------------------------------------------------------+
Output from code in #42 using GRASS 8.0.0RC2
> # after https://github.com/rsbivand/rgrass/issues/42#issuecomment-1040805090
> 
> library(rgrass7)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: (GRASS not running)
> my_GRASS <- "/usr/local/grass80"
> loc <- initGRASS(my_GRASS, home=tempdir(), gisDbase=".", 
+                  location="nc_basic_spm_grass7", mapset="user1", override=TRUE)
> 
> execGRASS("g.version")
GRASS 8.0.0RC2 (2022)
> execGRASS("g.gisenv")
GISDBASE='.';
LOCATION_NAME='nc_basic_spm_grass7';
MAPSET='user1';
GRASS_GUI='text';
> 
> execGRASS("g.region", raster="elevation")
> execGRASS("g.region", flags="p")
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228500
south:      215000
west:       630000
east:       645000
nsres:      10
ewres:      10
rows:       1350
cols:       1500
cells:      2025000
> 
> execGRASS("r.slope.aspect", flags="overwrite", elevation="elevation",
+           slope="slope", aspect="aspect")
 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete
> use_sp()
> u00 <- readRAST(c("elevation", "slope", "aspect")) # legacy
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Warning message:
In showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum NAD83 (High Accuracy Reference Network) in Proj4 definition
> object.size(u00) # in memory
48607248 bytes
> u0 <- rgrass7::read_RAST(c("elevation", "slope", "aspect"),
+                           return_format="SGDF") # re-written SGDF
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
Creating BIL support files...
Exporting raster as floating values (bytes=4)
 100%
> object.size(u0) # in memory
48607248 bytes
> library(sp)
> all.equal(u00[["elevation"]], u0[["elevation"]])
[1] TRUE
> u1 <- rgrass7:::read_RAST(c("elevation", "slope", "aspect"),
+                           return_format="terra") # re-written terra
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpMfAv7o/file82ea54a1122e.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpMfAv7o/file82ea2d74d6d1.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpMfAv7o/file82ea65297fac.grd> created.
> object.size(u1) # on file
1304 bytes
> library(terra)
terra 1.5.12
> u1
class       : SpatRaster 
dimensions  : 1350, 1500, 3  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_0=33.75 +lon_0=-79 +lat_1=36.1666666666667 +lat_2=34.3333333333333 +x_0=609601.22 +y_0=0 +ellps=GRS80 +units=m +no_defs 
sources     : file82ea54a1122e.grd  
              file82ea2d74d6d1.grd  
              file82ea65297fac.grd  
names       : elevation,     slope,    aspect 
min values  :  55.57879,   0.00000,   0.00000 
max values  : 156.32986,  38.68939, 360.00000 
> all.equal(u00[["elevation"]], c(values(u1[["elevation"]])))
[1] TRUE
> library(stars)
Loading required package: abind
Loading required package: sf
Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is TRUE
> u2 <- st_as_stars(u1) # coerce to stars
> u2
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                          Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
file82ea54a1122e.grd  55.57879 97.80552 118.1428 112.9484 130.4342 153.3789
dimension(s):
           from   to offset delta                       refsys point
x             1 1500 630000    10 +proj=lcc +lat_0=33.75 +l...    NA
y             1 1350 228500   -10 +proj=lcc +lat_0=33.75 +l...    NA
attributes    1    3     NA    NA                           NA    NA
                                    values x/y
x                                     NULL [x]
y                                     NULL [y]
attributes elevation, slope    , aspect       
> object.size(u2) # in memory
48611744 bytes
> all.equal(u00[["elevation"]], c(dplyr::pull(u2[,,,"elevation"])))
[1] TRUE
> u2a <- read_stars(u1@ptr$filenames, proxy=TRUE)
> u2a
stars_proxy object with 3 attributes in 3 file(s):
$file82ea54a1122e.grd
[1] "[...]/file82ea54a1122e.grd"

$file82ea2d74d6d1.grd
[1] "[...]/file82ea2d74d6d1.grd"

$file82ea65297fac.grd
[1] "[...]/file82ea65297fac.grd"

dimension(s):
  from   to offset delta                       refsys point values x/y
x    1 1500 630000    10 +proj=lcc +lat_0=33.75 +l...    NA   NULL [x]
y    1 1350 228500   -10 +proj=lcc +lat_0=33.75 +l...    NA   NULL [y]
> object.size(u2a) # on file
10328 bytes
> all.equal(u00[["elevation"]], c(st_as_stars(u2a)[[1]])) # moved to memory
[1] TRUE
> 
> 
> 
> # after https://github.com/rsbivand/rgrass/issues/42#issuecomment-1043311481
> 
> f <- system.file("ex/elev.tif", package="terra")
> r <- rast(f)
> plot(r, col=grDevices::terrain.colors(50))
> loc <- initGRASS(my_GRASS, home=tempdir(), SG=r, override=TRUE)
> 
> execGRASS("g.gisenv")
GISDBASE='/tmp/RtmpMfAv7o';
LOCATION_NAME='file82ea62f6e179';
MAPSET='file82ea49865a0b';
GRASS_GUI='text';
> 
> rgrass7:::write_RAST(r, "elev", flags="overwrite")
Importing raster map <elev>...
 100%
> execGRASS("r.info", map="elev")
 +----------------------------------------------------------------------------+
 | Map:      elev                           Date: Mon Feb 28 18:33:02 2022    |
 | Mapset:   file82ea49865a0b               Login of Creator: floris          |
 | Location: file82ea62f6e179                                                 |
 | DataBase: /tmp/RtmpMfAv7o                                                  |
 | Title:                                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    FCELL                Semantic label: (none)                |
 |   Rows:         90                                                         |
 |   Columns:      95                                                         |
 |   Total Cells:  8550                                                       |
 |        Projection: Latitude-Longitude                                      |
 |            N:  50:11:30N    S:  49:26:30N   Res: 0:00:30                   |
 |            E:      6:32E    W:   5:44:30E   Res: 0:00:30                   |
 |   Range of data:    min = 141  max = 547                                   |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal --overwrite input="/tmp/RtmpMfAv7o/file82ea3c994476.grd" o\   |
 |    utput="elev" memory=300 offset=0 num_digits=0                           |
 |                                                                            |
 +----------------------------------------------------------------------------+

> execGRASS("r.slope.aspect", flags="overwrite", elevation="elev", slope="slope", aspect="aspect")
 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete
> u1 <- rgrass7:::read_RAST(c("elev", "slope", "aspect"), return_format="terra")
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpMfAv7o/file82ea22025024.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpMfAv7o/file82ea7b75a041.grd> created.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
 100%
r.out.gdal complete. File </tmp/RtmpMfAv7o/file82ea2c69cdbd.grd> created.
> plot(u1[["elev"]], col=grDevices::terrain.colors(50))
> 
> 
> 
> # after https://github.com/rsbivand/rgrass/issues/42#issuecomment-1044769488
> 
> v1 <- terra::vect(system.file("ex/lux.shp", package = "terra"))
> v1
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 12, 6  (geometries, attributes)
 extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :  ID_1   NAME_1  ID_2   NAME_2  AREA       POP
 type        : <num>    <chr> <num>    <chr> <num>     <num>
 values      :     1 Diekirch     1 Clervaux   312 1.808e+04
                   1 Diekirch     2 Diekirch   218 3.254e+04
                   1 Diekirch     3  Redange   259 1.866e+04
> rgrass7:::write_VECT(v1, "lux")
Check if OGR layer <file82ea37d54a04> contains polygons...
 100%
Creating attribute table for layer <file82ea37d54a04>...
Default driver / database set to:
driver: sqlite
database: $GISDBASE/$LOCATION_NAME/$MAPSET/sqlite/sqlite.db
Importing 12 features (OGR layer <file82ea37d54a04>)...
 100%
-----------------------------------------------------
Registering primitives...

-----------------------------------------------------
Cleaning polygons
-----------------------------------------------------
Breaking polygons...
Breaking polygons (pass 1: select break points)...
 100%
Breaking polygons (pass 2: break at selected points)...
 100%
-----------------------------------------------------
Removing duplicates...
 100%
-----------------------------------------------------
Breaking boundaries...
 100%
-----------------------------------------------------
Removing duplicates...
 100%
-----------------------------------------------------
Cleaning boundaries at nodes...
 100%
-----------------------------------------------------
Merging boundaries...
 100%
-----------------------------------------------------
Changing boundary dangles to lines...
 100%
-----------------------------------------------------
Building areas...
-----------------------------------------------------
Changing boundary bridges to lines...
 100%
-----------------------------------------------------
Registering primitives...

Building areas...
 100%
Attaching islands...
 100%
-----------------------------------------------------
Finding centroids for OGR layer <file82ea37d54a04>...
 100%
-----------------------------------------------------
Writing centroids...
 100%
-----------------------------------------------------
12 input polygons
Total area: 2.56481E+09 (12 areas)
-----------------------------------------------------
Copying features...
 100%
Building topology for vector map <lux@file82ea49865a0b>...
Registering primitives...

Building areas...
 100%
Attaching islands...
 100%
Attaching centroids...
 100%
Warning message:
In x@ptr$write(filename, layer, filetype, overwrite[1], options) :
  GDAL Message 6: dataset /tmp/RtmpMfAv7o/file82ea37d54a04.gpkg does not support layer creation option ENCODING
> execGRASS("v.info", map="lux")
 +----------------------------------------------------------------------------+
 | Name:            lux                                                       |
 | Mapset:          file82ea49865a0b                                          |
 | Location:        file82ea62f6e179                                          |
 | Database:        /tmp/RtmpMfAv7o                                           |
 | Title:                                                                     |
 | Map scale:       1:1                                                       |
 | Name of creator: floris                                                    |
 | Organization:                                                              |
 | Source date:     Mon Feb 28 18:33:03 2022                                  |
 | Timestamp (first layer): none                                              |
 |----------------------------------------------------------------------------|
 | Map format:      native                                                    |
 |----------------------------------------------------------------------------|
 |   Type of map: vector (level: 2)                                           |
 |                                                                            |
 |   Number of points:       0               Number of centroids:  12         |
 |   Number of lines:        0               Number of boundaries: 33         |
 |   Number of areas:        12              Number of islands:    1          |
 |                                                                            |
 |   Map is 3D:              No                                               |
 |   Number of dblinks:      1                                                |
 |                                                                            |
 |   Projection: Latitude-Longitude                                           |
 |                                                                            |
 |               N:   50:10:53.83758N    S:  49:26:52.106316N                 |
 |               E:   6:31:41.707632E    W:    5:44:38.90454E                 |
 |                                                                            |
 |   Digitization threshold: 0                                                |
 |   Comment:                                                                 |
 |                                                                            |
 +----------------------------------------------------------------------------+

man/initGRASS.Rd Outdated Show resolved Hide resolved
R/rast_link.R Outdated Show resolved Hide resolved
Copy link
Collaborator

@florisvdh florisvdh left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks fine; nice additions! No problems reported by R CMD check using this rgrass7 under latest releases GRASS 7.8.7 and 8.0.1, in Linux Mint 20.

For the examples in initGRASS.Rd, note they will only run (and be checked) provided that the GRASS_INSTALLATION variable is set. Although having the possibility to externally set the GRASS installation directory is an interesting addition IMHO, you may still want to add the following in initGRASS.Rd, for the usual case where GRASS_INSTALLATION is empty (so that the code would be run during a regular check with the nc_basic_spm_grass7 location):

rgrass/man/readRAST.Rd

Lines 60 to 62 in 5bfe117

run <- FALSE
if (nchar(Sys.getenv("GISRC")) > 0 &&
read.dcf(Sys.getenv("GISRC"))[1,"LOCATION_NAME"] == "nc_basic_spm_grass7") run <- TRUE

Further, I've sent PR #46 that goes to this rastNG branch; if accepted that one is to be merged first in order to add it to current PR. (Note, I normally always commit to a separate branch, in order to avoid conflicting pushes to the same branch by different contributors).

@rsbivand
Copy link
Owner Author

rsbivand commented Mar 3, 2022

@florisvdh I don't think that enabling the initGRASS() examples when running R in an active GRASS environment is sensible, because GRASS is already running. There is no way to test initGRASS() (GRASS in R) and regular (R in GRASS) at the same time, I think.

readVECT.Rd: remove layers from GRASS mapset
@florisvdh
Copy link
Collaborator

There is no way to test initGRASS() (GRASS in R) and regular (R in GRASS) at the same time, I think.

Could be, depending on the circumstances, although I didn't run into problems so far. Rerunning initGRASS() in a R-in-GRASS session just resets the needed env variables and seems to work well, using the same GRASS version, although I wouldn't try/recommend to call another GRASS version that way. The examples at the bottom of #45 (review) were run as R within GRASS and also contain a change of the location. Rerunning some lines, from R within GRASS 8.0.1 (actually RStudio, see GRASS session code):

Invoking RStudio from GRASS (click to expand)
$ /usr/local/grass8.0.1/bin/grass 
Starting GRASS GIS...
Cleaning up temporary files...

          __________  ___   __________    _______________
         / ____/ __ \/   | / ___/ ___/   / ____/  _/ ___/
        / / __/ /_/ / /| | \__ \\_  \   / / __ / / \__ \
       / /_/ / _, _/ ___ |___/ /__/ /  / /_/ // / ___/ /
       \____/_/ |_/_/  |_/____/____/   \____/___//____/

Welcome to GRASS GIS 8.0.1
GRASS GIS homepage:                      https://grass.osgeo.org
This version running through:            Bash Shell (/bin/bash)
Help is available with the command:      g.manual -i
See the licence terms with:              g.version -c
See citation options with:               g.version -x
Start the GUI with:                      g.gui wxpython
When ready to quit enter:                exit

To run a command as administrator (user "root"), use "sudo <command>".
See "man sudo_root" for details.

GRASS nc_basic_spm_grass7/PERMANENT:git_repositories > rstudio 2> /dev/null &
[1] 6947
GRASS nc_basic_spm_grass7/PERMANENT:git_repositories > 
R code and output (click to expand)
> library(rgrass7)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: GRASS 8.0.1 (2022)
and location: nc_basic_spm_grass7
> execGRASS("g.version")
GRASS 8.0.1 (2022)
> execGRASS("g.gisenv")
GISDBASE='/home/floris/grassdata';
LOCATION_NAME='nc_basic_spm_grass7';
MAPSET='PERMANENT';
GUI='text';
PID='6927';
>
>
> setwd("~/ContinuLegen/rgrass_test/")
> getwd()
[1] "/media/floris/DATA/ContinuLegen/rgrass_test"
> my_GRASS <- "/usr/local/grass8.0.1/grass80"
> loc <- initGRASS(my_GRASS, home=tempdir(), gisDbase=".", 
+                  location="nc_basic_spm_grass7", mapset="user1", override=TRUE)
> execGRASS("g.version")
GRASS 8.0.1 (2022)
> execGRASS("g.gisenv")
GISDBASE='.';
LOCATION_NAME='nc_basic_spm_grass7';
MAPSET='user1';
GRASS_GUI='text';
>
>
> library(terra)
terra 1.5.12
> f <- system.file("ex/elev.tif", package="terra")
> r <- rast(f)
> loc <- initGRASS(my_GRASS, home=tempdir(), SG=r, override=TRUE)
> execGRASS("g.gisenv")
GISDBASE='/tmp/grass8-floris-6922/RtmpvZDi8a';
LOCATION_NAME='file1b615d871075';
MAPSET='file1b616c223a0d';
GRASS_GUI='text';

@rsbivand
Copy link
Owner Author

rsbivand commented Mar 3, 2022

I think I'd prefer to write a short vignette rather than use examples to demonstrate the different wayss of using the package, as a new PR. I just marked the old functions as deprecated, which may alert users to forthcoming changes. Are we OK to merge?

@florisvdh
Copy link
Collaborator

Yes, a vignette is the best way.

I just marked the old functions as deprecated

Not pushed that commit yet?

Are we OK to merge?

I didn't encounter problems and totally trust you!

@rsbivand
Copy link
Owner Author

rsbivand commented Mar 4, 2022

Thanks, now pushed! And that's why trusting totally implies trusting collaborative processes!

R/bin_link.R Outdated Show resolved Hide resolved
rsbivand and others added 2 commits March 4, 2022 11:37
Co-authored-by: Floris Vanderhaeghe <floris.vanderhaeghe@inbo.be>
@rsbivand
Copy link
Owner Author

rsbivand commented Mar 4, 2022

Thanks, I updated the other deprecation messages.

@florisvdh
Copy link
Collaborator

Thanks, all looks good!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants