Re: [R-sig-Geo] v.split.length (GRASS) in R

2016-11-23 Thread julian magezi via R-sig-Geo
Hello members,
Which package in R does nonparametric statistical testing?
Thanks
 Julian Ijumulana, MSc (Geoinformatics)Land surveyor, GIS/RS Specialist & 
Assistant Lecturer, Department of Transportation and Geotechnical 
Engineering,College of Engineering and Technology,University of Dar es 
Salaam,Box 35131,Dar es SalaamTANZANIA
Mob: +255 767675416

 

On Wednesday, 23 November 2016, 20:49, Manuel Spínola 
 wrote:
 

 Thank you very much Florian.

Do you know any R function/package to do the same without relying on GRASS?

Manuel

2016-11-23 7:28 GMT-06:00 Florian Betz :

> Assuming you are using GRASS 7, you can use the rgrass7 package. First
> thing to do is to initialize the R-GRASS connection, then you can run the
> module using execGRASS.
>
> library(rgrass7)
> #Arguments depending on your system
> initGRASS(gisBase ="C:/Program Files/GRASS GIS 7.0.5",home=tempdir(),
>          gisDbase="Path to the folder with your location, location="Your
> GRASS location",mapset="Your GRASS mapset",override=TRUE)
> #See the help of v.split.lenght for the necessary function arguments
> execGRASS("v.split.length", argument1=, argument2=, ...)
>
> Regards,
>
> Flo
>
>
> Am 23.11.2016 um 14:15 schrieb Manuel Spínola:
>
>> Dear list members,
>>
>> Is it possible to run v.split.length from GRASS in R? Or is there anyway
>> to
>> split a SpatialLinesDataFrame to shorter segments by length in R?
>>
>> Best,
>>
>> Manuel
>>
>>
>


-- 
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspin...@una.cr 
mspinol...@gmail.com
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río 
Institutional website: ICOMVIS 

    [[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

   
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] v.split.length (GRASS) in R

2016-11-23 Thread Manuel Spínola
Thank you very much Jannes.

Does the example works with QGIS 2.18 or should I do some setting in QGIS
before running the code?

Manuel

2016-11-23 8:45 GMT-06:00 "Jannes Münchow" :

> To achieve this, you can also use the RQGIS-package in conjunction with
> GRASS given you have installed QGIS along with GRASS (for more information
> have a look at vignette("install_guide", package = "RQGIS"):
>
> # construct a SpatialPointsDataFrame (which will be accepted as input by
> # run_qgis)
> library("sp")
> # from the sp vignette:
> l1 <- cbind(c(1, 2, 3), c(3, 2, 2))
> rownames(l1) <- letters[1:3]
> l1a <- cbind(l1[, 1] + 0.05, l1[, 2] + 0.05)
> rownames(l1a) <- letters[1:3]
> l2 <- cbind(c(1, 2, 3), c(1, 1.5, 1))
> rownames(l2) <- letters[1:3]
> Sl1 <- Line(l1)
> Sl1a <- Line(l1a)
> Sl2 <- Line(l2)
> S1 <- Lines(list(Sl1, Sl1a), ID = "a")
> S2 <- Lines(list(Sl2), ID = "b")
> Sl <- SpatialLines(list(S1, S2))
> # convert it to a SpatialLinesDataFrame
> Sl <- SpatialLinesDataFrame(Sl, data = data.frame(1:2), match.ID = FALSE)
> proj4string(Sl) <- CRS("+proj=longlat +datum=WGS84")
>
> # Now use RQGIS
> library("RQGIS")
> # indicate where QGIS is installed on your computer
> qgis_env <- set_env("C:/OSGeo4W64/")
> args <- get_args_man("grass7:v.split.length", qgis_env = qgis_env,
>  options = TRUE)
> # have a look at the GRASS online help
> open_help("grass7:v.split.length", qgis_env = qgis_env)
>
> # specify the necessary arguments
> args$input <- Sl
> # here length corresponds to one decimal degree
> args$length <- "1"
> args$output <- file.path(tempdir(), "out.shp")
> # load the output directly into R again
> out <- run_qgis(alg = "grass7:v.split.length", params = args,
> load_output = args$output,
> qgis_env = qgis_env)
> length(Sl)  # 2 line objects
> length(out)  #  9 line objects
> # Have a look at the output
> plot(out, col = rep(c("blue", "green", "black"), 3))
>
> Cheers,
>
> Jannes
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo




-- 
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspin...@una.cr 
mspinol...@gmail.com
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río 
Institutional website: ICOMVIS 

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] v.split.length (GRASS) in R

2016-11-23 Thread Manuel Spínola
Thank you very much Florian.

Do you know any R function/package to do the same without relying on GRASS?

Manuel

2016-11-23 7:28 GMT-06:00 Florian Betz :

> Assuming you are using GRASS 7, you can use the rgrass7 package. First
> thing to do is to initialize the R-GRASS connection, then you can run the
> module using execGRASS.
>
> library(rgrass7)
> #Arguments depending on your system
> initGRASS(gisBase ="C:/Program Files/GRASS GIS 7.0.5",home=tempdir(),
>   gisDbase="Path to the folder with your location, location="Your
> GRASS location",mapset="Your GRASS mapset",override=TRUE)
> #See the help of v.split.lenght for the necessary function arguments
> execGRASS("v.split.length", argument1=, argument2=, ...)
>
> Regards,
>
> Flo
>
>
> Am 23.11.2016 um 14:15 schrieb Manuel Spínola:
>
>> Dear list members,
>>
>> Is it possible to run v.split.length from GRASS in R? Or is there anyway
>> to
>> split a SpatialLinesDataFrame to shorter segments by length in R?
>>
>> Best,
>>
>> Manuel
>>
>>
>


-- 
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspin...@una.cr 
mspinol...@gmail.com
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río 
Institutional website: ICOMVIS 

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster predict with gnls model

2016-11-23 Thread Johannes Schumacher via R-sig-Geo
Dear List,

I have trouble with the predict function (raster package) when using a gnls 
model (nonlinear model fit using generalized least squares, from the nlme 
package).


I'm using R version 3.3.1 (2016-06-21), Platform: x86_64-w64-mingw32/x64 
(64-bit), and raster package version 2.5-8.
I want to make predictions based on three continuous variables and one 
categorical variable. The predictors are layers in a RasterStack object (see 
description below: s20, s20_excludefactor). Using other models than gnls works. 
I tested three models (my_nls, my_lm, my_gnls, see description below):

1) Using nonlinear least squares (nls) fit excluding the factor as predictor 
works fine:
v20 <- predict(object=s20_excludefactor, model=my_nls) 
2) Using linear model (lm) including the factor as predictor works as well, 
except the warning:
“not sure if the correct factor levels are used here”. 
v20 <- predict(object=s20, model=my_lm) #  or
v20 <- predict(object=s20, model=my_lm, na.rm=T, 
factors=list(Gruppe=c("1","2","3")))

3) Using the glns model produces error:
“Error in p %*% beta[pmap[[nm]]] : non-conformable arguments“
v20 <- predict(object=s20, model=my_gnls) #  or
v20 <- predict(object=s20, model=my_gnls, na.rm=T, 
factors=list(Gruppe=c("1","2","3")))


I don't know what to do about the error in 3). Any help is highly appreciated! 
Best regards,
Johannes

> s20_excludefactor (used in my_nls)
class   : RasterStack 
dimensions  : 50, 50, 2500, 3  (nrow, ncol, ncell, nlayers)
resolution  : 20, 20  (x, y)
extent  : 339, 3391000, 5296000, 5297000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 
+datum=potsdam +units=m +no_defs +ellps=bessel 
+towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 
names   : x1, x2, x3 
min values  :  0.,  0.,207.3915 
max values  : 243534.0084,  1.,220.1788 


> s20 (used in my_lm and my_gnls)
class   : RasterStack 
dimensions  : 50, 50, 2500, 4  (nrow, ncol, ncell, nlayers)
resolution  : 20, 20  (x, y)
extent  : 339, 3391000, 5296000, 5297000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 
+datum=potsdam +units=m +no_defs +ellps=bessel 
+towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 
names   : x1, x2, x3,  Gruppe 
min values  :  0.,  0.,207.3915,  1. 
max values  : 243534.0084,  1.,220.1788,  3. 


> my_nls
Nonlinear regression model
model: y ~ a/(1 + exp(b0 + b1 * x1 + b2 * x2 + b3 * x3))
data: datapply
a b0 b1 b2 b3 
1.500e+03  4.544e+00 -9.134e-06 -1.084e+00 -6.121e-04 
weighted residual sum-of-squares: 56.55


> my_lm
Call:
lm(formula = y ~ x1 + x2 + x3 + Gruppe, data = datapply)
Coefficients:
(Intercept) x1 x2 x3  Gruppe2  Gruppe3 
8.208e+012.641e-03   -3.463e+028.357e-026.749e+018.759e+01 


> my_gnls
Generalized nonlinear least squares fit
Model: y ~ a/(1 + exp(b0 + b1 * x1 + b2 * x2 + b3 * x3)) 
Data: datapply 
Log-likelihood: -4305.82
Coefficients:
a b0 b1 b2 b3.(Intercept) b3.Gruppe2
 b3.Gruppe3 
1.312412e+03   4.410668e+00  -1.001456e-05  -1.082276e+00  -7.200440e-05  
-4.635569e-04  -5.423707e-04 
Variance function:
Structure: Power of variance covariate
Formula: ~x1 
Parameter estimates:
power 
0.9695105 

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] R - QGIS integration

2016-11-23 Thread MacQueen, Don
And if Kent is correct, then there is also:


subplot  package:Hmisc  R Documentation

Embed a new plot within an existing plot

Description:

 Subplot will embed a new plot within an existing plot at the
 coordinates specified (in user units of the existing plot).



(Which would not require learning ggplot2)
-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 11/22/16, 6:05 AM, "R-sig-Geo on behalf of Kent Johnson"
 wrote:

>If I understand correctly your question is not about QGIS but about using
>R
>to make maps with inset charts. There are several examples of charts on
>maps here:
>http://stackoverflow.com/questions/10368180/plotting-pie-graphs-on-map-in-
>ggplot
>
>Also ggmap + ggmap::inset looks like it will do what you want, see the
>ggplot2::annotation_custom example for a little help.
>
>Kent
>
>On Mon, Nov 21, 2016 at 6:00 AM,  wrote:
>
>
>> -
>>
>> Message: 7
>> Date: Mon, 21 Nov 2016 10:40:26 +0100
>> From: Domagoj Culinovic 
>> To: R-sig-Geo@r-project.org
>> Subject: [R-sig-Geo] R - QGIS integration
>> Message-ID:
>> > gmail.com>
>> Content-Type: text/plain; charset="UTF-8"
>>
>> Hi,
>> I have Polygon and point data in SHP format which i have use in QGIS.
>>also
>> i have several csv files where i have different statistical data about
>> points and polygons in SHP files.
>> This csv data are generated every week.
>> How to make R script which can make georeferenced charts, diagrams or
>> something like heathmap i R, and return picture or charts at theit
>>position
>> on the map?
>>
>> Regards,
>> Domagoj
>>
>> [[alternative HTML version deleted]]
>>
>>
>>
>> --
>>
>> Subject: Digest Footer
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> --
>>
>> End of R-sig-Geo Digest, Vol 159, Issue 15
>> **
>>
>
>   [[alternative HTML version deleted]]
>
>___
>R-sig-Geo mailing list
>R-sig-Geo@r-project.org
>https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


Re: [R-sig-Geo] Calculate the area of pixel in square kilometers from a raster in lat/long projection

2016-11-23 Thread Isaque Daniel
Hi Thiago,


The function calculate the area by pixel "like raster::area", but, considering 
the distortions of projection and transforming from degree to hectares (or 
square kilometers, if someone need).

The problem is: how to compute the area for raster maps during the interactive 
process and compare as the values used to define the classification (chosing 
the thresholds), because the raster::area return the area in the projection of 
raster, and for huge datasets or multiple aplications in sequence (or both 
togheter, as in my case), is computational hard to reproject raster ou poor to 
work with "mean value of pixel" in large rasters..

So, I develop this simple function.
Best
Isaque


--
Eng. Agr. Isaque Daniel Rocha Eberhardt
Mestre em Sensoriamento Remoto - Instituto Nacional de Pesquisas Espaciais 
(INPE)
Doutorando em Transportes - Universidade de Bras??lia (UNB)
Mobile: +55 (061) 99015658
--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Bras??lia University (UNB)



De: Thiago V. dos Santos 
Enviado: quarta-feira, 23 de novembro de 2016 15:56
Para: Isaque Daniel; r-sig-geo
Assunto: Re: [R-sig-Geo] Calculate the area of pixel in square kilometers from 
a raster in lat/long projection

raster::area would calculate the area in square kilometers, and from there one 
can convert to hectare or square meters.

I was wondering: how do those calculations differ from raster::area's?
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota


On Tuesday, November 22, 2016 8:15 PM, Isaque Daniel 
 wrote:



Hi there,
I found a example in stackoverflow of distance calculation and convertion, so I 
modificate that to compute the area by pixel in hectares using the raster in 
degrees.

#
# Function to compute the area  in hectares

measureAreahec <- function(imageRaster) {
  R <- 6378.137# radius of earth in Km
  areaImagetemp <- area(ceiImg)
  xresVal <- xres(areaImagetemp)
  yresVal <- yres(areaImagetemp)
  coords <- data.frame(coordinates(areaImagetemp))
  colnames(coords) <- c("lon","lat")
  lon1 <- coords$lon - (xresVal/2)
  lon2 <- coords$lon + (xresVal/2)
  lat1 <- coords$lat - (yresVal/2)
  lat2 <- coords$lat + (yresVal/2)

  # width
  dLat <- (lat1-lat1)*pi/180
  dLon <- (lon2-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlongitude <- R * c * 1000

  # heigth
  dLat <- (lat1-lat2)*pi/180
  dLon <- (lon1-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlatitude <- R * c * 1000
  return (dlatitude * dlongitude / 1)# distance 
in meters
}
##


It's a very simple one, but ,maybe can help somebody.

Best
Isaque


--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Brazlia University (UNB)



De: R-sig-Geo  em nome de Isaque Daniel 

Enviado: ter�-feira, 22 de novembro de 2016 17:48
Para: r-sig-geo
Assunto: [R-sig-Geo] Calculate the area of pixel in square kilometers from a 
raster in lat/long projection

Hi,


I'm looking for some advice to convert the area calculated by raster::area() in 
a raster image as the projection in lat/long.

[[elided Hotmail spam]]
Best
Isaque


--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Brazilia University (UNB)

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
R-sig-Geo Info Page
stat.ethz.ch
R-sig-Geo -- R Special Interest Group on using Geographical data and Mapping 
About R-sig-Geo


R-sig-Geo Info Page
R-sig-Geo Info Page
stat.ethz.ch
R-sig-Geo -- R Special Interest Group on using Geographical data and Mapping 
About R-sig-Geo


stat.ethz.ch
R-sig-Geo -- R Special Interest Group on using Geographical data and Mapping 
About 

Re: [R-sig-Geo] Calculate the area of pixel in square kilometers from a raster in lat/long projection

2016-11-23 Thread Thiago V. dos Santos via R-sig-Geo
raster::area would calculate the area in square kilometers, and from there one 
can convert to hectare or square meters.

I was wondering: how do those calculations differ from raster::area's? 
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota


On Tuesday, November 22, 2016 8:15 PM, Isaque Daniel 
 wrote:



Hi there,
I found a example in stackoverflow of distance calculation and convertion, so I 
modificate that to compute the area by pixel in hectares using the raster in 
degrees.

#
# Function to compute the area  in hectares

measureAreahec <- function(imageRaster) {
  R <- 6378.137# radius of earth in Km
  areaImagetemp <- area(ceiImg)
  xresVal <- xres(areaImagetemp)
  yresVal <- yres(areaImagetemp)
  coords <- data.frame(coordinates(areaImagetemp))
  colnames(coords) <- c("lon","lat")
  lon1 <- coords$lon - (xresVal/2)
  lon2 <- coords$lon + (xresVal/2)
  lat1 <- coords$lat - (yresVal/2)
  lat2 <- coords$lat + (yresVal/2)

  # width
  dLat <- (lat1-lat1)*pi/180
  dLon <- (lon2-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlongitude <- R * c * 1000

  # heigth
  dLat <- (lat1-lat2)*pi/180
  dLon <- (lon1-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlatitude <- R * c * 1000
  return (dlatitude * dlongitude / 1)# distance 
in meters
}
##


It's a very simple one, but ,maybe can help somebody.

Best
Isaque


--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Brazlia University (UNB)



De: R-sig-Geo  em nome de Isaque Daniel 

Enviado: ter�-feira, 22 de novembro de 2016 17:48
Para: r-sig-geo
Assunto: [R-sig-Geo] Calculate the area of pixel in square kilometers from a 
raster in lat/long projection

Hi,


I'm looking for some advice to convert the area calculated by raster::area() in 
a raster image as the projection in lat/long.

[[elided Hotmail spam]]
Best
Isaque


--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Brazilia University (UNB)

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
R-sig-Geo Info Page
stat.ethz.ch
R-sig-Geo -- R Special Interest Group on using Geographical data and Mapping 
About R-sig-Geo





[[alternative HTML version deleted]]


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] v.split.length (GRASS) in R

2016-11-23 Thread Jannes Münchow
To achieve this, you can also use the RQGIS-package in conjunction with GRASS 
given you have installed QGIS along with GRASS (for more information have a 
look at vignette("install_guide", package = "RQGIS"):

# construct a SpatialPointsDataFrame (which will be accepted as input by
# run_qgis)
library("sp")
# from the sp vignette:
l1 <- cbind(c(1, 2, 3), c(3, 2, 2))
rownames(l1) <- letters[1:3]
l1a <- cbind(l1[, 1] + 0.05, l1[, 2] + 0.05)
rownames(l1a) <- letters[1:3]
l2 <- cbind(c(1, 2, 3), c(1, 1.5, 1))
rownames(l2) <- letters[1:3]
Sl1 <- Line(l1)
Sl1a <- Line(l1a)
Sl2 <- Line(l2)
S1 <- Lines(list(Sl1, Sl1a), ID = "a")
S2 <- Lines(list(Sl2), ID = "b")
Sl <- SpatialLines(list(S1, S2))
# convert it to a SpatialLinesDataFrame
Sl <- SpatialLinesDataFrame(Sl, data = data.frame(1:2), match.ID = FALSE)
proj4string(Sl) <- CRS("+proj=longlat +datum=WGS84")

# Now use RQGIS
library("RQGIS")
# indicate where QGIS is installed on your computer
qgis_env <- set_env("C:/OSGeo4W64/")
args <- get_args_man("grass7:v.split.length", qgis_env = qgis_env, 
 options = TRUE)
# have a look at the GRASS online help
open_help("grass7:v.split.length", qgis_env = qgis_env)

# specify the necessary arguments
args$input <- Sl
# here length corresponds to one decimal degree
args$length <- "1"
args$output <- file.path(tempdir(), "out.shp")
# load the output directly into R again
out <- run_qgis(alg = "grass7:v.split.length", params = args,
load_output = args$output,
qgis_env = qgis_env)
length(Sl)  # 2 line objects
length(out)  #  9 line objects
# Have a look at the output
plot(out, col = rep(c("blue", "green", "black"), 3))
 
Cheers, 

Jannes

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] rgrass7 : Error in parseGRASS

2016-11-23 Thread Michael DELORME
Thanks for sharing ; I'll look into it.
Cordially

Le 23/11/2016 09:58, Chris Reudenbach a écrit :
> Michael,
>
>
> using GRASS from R is a bit cumbersome. GRASS (and this is for sure
> one of it strength ) is extremely rigid with projections etc. Due to
> this you should setup a correct structure. I use a setup that is
> similar to the below snippet.
>
> cheers Chris
>
>
> library(rgrass7)
> library(sp)
>
> # setup a temp workingdir
> working.dir<- "~/tmp/"
>
> # get meuse data
>
> data(meuse)
> data(meuse.grid)
>
> # georeference the meuse grid data
> coordinates(meuse.grid) <- ~x+y
> proj4string(meuse.grid) <- CRS("+init=epsg:28992")
> gridded(meuse.grid) <- TRUE
>
> # get a first cellsize/pixel resolution for GRASS
> resolution <- meuse.grid@grid@cellsize[1]
>
> # georeference the meuse data
> coordinates(meuse) <- ~x+y
> proj4string(meuse) <- as.character(CRS("+init=epsg:28992"))
>
> # get projection, proj4 string and extent for GRASS
> projection<-(strsplit(meuse@proj4string@projargs,split = " "))
> proj4<- paste(projection[[1]][2:length(unlist(projection))], collapse
> = ' ')
> xmax<-meuse@bbox[3]
> xmin<-meuse@bbox[1]
> ymax<-meuse@bbox[4]
> ymin<-meuse@@bbox[2]
>
> # create and set working directory
> if (!file.exists(file.path(working.dir,"run"))){
>   dir.create(file.path(working.dir,"run"),recursive = TRUE)
> }
> setwd(file.path( working.dir,"run"))
>
>
>  set up GRASS RUNTIME environment
>
> # define the GRASS executable path
> if(Sys.info()["sysname"] == "Windows"){
>   grass.gis.base<-'C:/OSGeo4W64/apps/grass/grass-7.0.5'
> }else {
>   grass.gis.base<-'/usr/lib/grass72'
> }
>
> # set path for optional GRASS addons
> Sys.setenv(GRASS_ADDON_PATH="~/.grass7/addons")
>
> # create the TEMPORARY GRASS location
> rgrass7::initGRASS(gisBase=grass.gis.base,
>home=tempdir(),
>mapset='PERMANENT',
>override=TRUE
> )
>
> # assign GRASS projection according to data set
> rgrass7::execGRASS('g.proj',
>flags=c('c','quiet'),
>proj4=proj4
> )
>
> # assign GRASS extent and resolution
> rgrass7::execGRASS('g.region',
>flags=c('quiet'),
>n=as.character(ymax),
>s=as.character(ymin),
>e=as.character(xmax),
>w=as.character(xmin),
>res=as.character(resolution)
> )
>
>
> #   now do GRASS STUFF
>
> writeVECT(meuse,"meuse")
> execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 double
> precision")
>
>
> On 23.11.2016 07:38, Michael DELORME wrote:
>> Thanks for your insight.
>> I'll use a stand-alone GRASS install.
>>
>> Cordially
>>
>> Le 22/11/2016 15:13, Roger Bivand a écrit :
>>> On Tue, 22 Nov 2016, Michael DELORME wrote:
>>>
 Thanks for investigating.

 I'll try a standalone version of GRASS if needed.

 Here is the traceback

> traceback()
 4: stop(paste(cmd, "not parsed"))
 3: parseGRASS(cmd, legacyExec = legacyExec)
 2: doGRASS(cmd, flags = flags, ..., parameters = parameters, echoCmd =
 echoCmd, legacyExec = legacyExec)
 1: execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2
 double
 precision")

 My C:\OSGeo4W64\apps\grass\grass-7.0.5\bin lists *.exe or *.bat and
 indeed the *.bat command are those which don't work (r.mask.bat,
 v.db.addcolumn.bat,...). The batch calls Python scripts ; for example
 v.db.addcolumn.bat is :

 @"%GRASS_PYTHON%" "%GISBASE%/scripts/v.db.addcolumn.py" %*

 So I guess something must be wrong in my Python installation... (it
 works however from within a GRASS shell)
>>> No, the logic for Windows and initGRASS has only ever been tested with
>>> stand-alone GRASS. I would guess that in OSGeo4W the path variables
>>> are wrong. I cannot even find out how to run initGRASS in R but
>>> outside the OSGeo4W shell (doesn't find iconv.dll) or inside the
>>> OSGeo4W shell (cannot create a temporary GRASS location). Unless you
>>> or others really need OSGeo4W, I suggest asking others to implement
>>> that - very bulky and fragile.
>>>
>>> Roger
>>>
 Cordially


 Le 22/11/2016 14:30, Roger Bivand a écrit :
> On Tue, 22 Nov 2016, Michael DELORME wrote:
>
>> Thanks for replying
>>
>> Here is a reproducible example (change your GRASS directory of
>> course) :
>>
>> library(rgrass7)
>> library(sp)
>> initGRASS("C:/OSGeo4W64/apps/grass/grass-7.0.5", home = tempdir(),
>> override = TRUE)
>> data(meuse)
>> coordinates(meuse) <- ~x+y
>> writeVECT(meuse, "meuse")
>> execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 double
>> precision")
> This appears to run correctly on Windows 7 with GRASS 7.0.5 Windows
> stand-alone. Do you need to use OSGeo4W - the internal logic used for
> finding out where the different versions keep their executables and

Re: [R-sig-Geo] rgrass7 : Error in parseGRASS

2016-11-23 Thread Roger Bivand

Chris,

Thanks for this. Could you please explain how you use R in the OSGeo4W64 
context? Where are you starting R and how does it find the OSGeo4W 
components that are not on default PATH settings? In order to correct 
rgrass7, I need to know how to start R under OSGeo4W64 outside GRASS.


Roger

On Wed, 23 Nov 2016, Chris Reudenbach wrote:


Michael,


using GRASS from R is a bit cumbersome. GRASS (and this is for sure one of it 
strength ) is extremely rigid with projections etc. Due to this you should 
setup a correct structure. I use a setup that is similar to the below 
snippet.


cheers Chris


library(rgrass7)
library(sp)

# setup a temp workingdir
working.dir<- "~/tmp/"

# get meuse data

data(meuse)
data(meuse.grid)

# georeference the meuse grid data
coordinates(meuse.grid) <- ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) <- TRUE

# get a first cellsize/pixel resolution for GRASS
resolution <- meuse.grid@grid@cellsize[1]

# georeference the meuse data
coordinates(meuse) <- ~x+y
proj4string(meuse) <- as.character(CRS("+init=epsg:28992"))

# get projection, proj4 string and extent for GRASS
projection<-(strsplit(meuse@proj4string@projargs,split = " "))
proj4<- paste(projection[[1]][2:length(unlist(projection))], collapse = ' ')
xmax<-meuse@bbox[3]
xmin<-meuse@bbox[1]
ymax<-meuse@bbox[4]
ymin<-meuse@@bbox[2]

# create and set working directory
if (!file.exists(file.path(working.dir,"run"))){
 dir.create(file.path(working.dir,"run"),recursive = TRUE)
}
setwd(file.path( working.dir,"run"))


 set up GRASS RUNTIME environment

# define the GRASS executable path
if(Sys.info()["sysname"] == "Windows"){
 grass.gis.base<-'C:/OSGeo4W64/apps/grass/grass-7.0.5'
}else {
 grass.gis.base<-'/usr/lib/grass72'
}

# set path for optional GRASS addons
Sys.setenv(GRASS_ADDON_PATH="~/.grass7/addons")

# create the TEMPORARY GRASS location
rgrass7::initGRASS(gisBase=grass.gis.base,
   home=tempdir(),
   mapset='PERMANENT',
   override=TRUE
)

# assign GRASS projection according to data set
rgrass7::execGRASS('g.proj',
   flags=c('c','quiet'),
   proj4=proj4
)

# assign GRASS extent and resolution
rgrass7::execGRASS('g.region',
   flags=c('quiet'),
   n=as.character(ymax),
   s=as.character(ymin),
   e=as.character(xmax),
   w=as.character(xmin),
   res=as.character(resolution)
)


#   now do GRASS STUFF

writeVECT(meuse,"meuse")
execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 double
precision")


On 23.11.2016 07:38, Michael DELORME wrote:

 Thanks for your insight.
 I'll use a stand-alone GRASS install.

 Cordially

 Le 22/11/2016 15:13, Roger Bivand a écrit :
>  On Tue, 22 Nov 2016, Michael DELORME wrote:
> 
> >  Thanks for investigating.
> > 
> >  I'll try a standalone version of GRASS if needed.
> > 
> >  Here is the traceback
> > 
> > >  traceback()

> >  4: stop(paste(cmd, "not parsed"))
> >  3: parseGRASS(cmd, legacyExec = legacyExec)
> >  2: doGRASS(cmd, flags = flags, ..., parameters = parameters, echoCmd =
> >  echoCmd, legacyExec = legacyExec)
> >  1: execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 
> >  double

> >  precision")
> > 
> >  My C:\OSGeo4W64\apps\grass\grass-7.0.5\bin lists *.exe or *.bat and

> >  indeed the *.bat command are those which don't work (r.mask.bat,
> >  v.db.addcolumn.bat,...). The batch calls Python scripts ; for example
> >  v.db.addcolumn.bat is :
> > 
> >  @"%GRASS_PYTHON%" "%GISBASE%/scripts/v.db.addcolumn.py" %*
> > 
> >  So I guess something must be wrong in my Python installation... (it

> >  works however from within a GRASS shell)
>  No, the logic for Windows and initGRASS has only ever been tested with
>  stand-alone GRASS. I would guess that in OSGeo4W the path variables
>  are wrong. I cannot even find out how to run initGRASS in R but
>  outside the OSGeo4W shell (doesn't find iconv.dll) or inside the
>  OSGeo4W shell (cannot create a temporary GRASS location). Unless you
>  or others really need OSGeo4W, I suggest asking others to implement
>  that - very bulky and fragile.
> 
>  Roger
> 
> >  Cordially
> > 
> > 
> >  Le 22/11/2016 14:30, Roger Bivand a écrit :

> > >  On Tue, 22 Nov 2016, Michael DELORME wrote:
> > > 
> > > >  Thanks for replying
> > > > 
> > > >  Here is a reproducible example (change your GRASS directory of

> > > >  course) :
> > > > 
> > > >  library(rgrass7)

> > > >  library(sp)
> > > >  initGRASS("C:/OSGeo4W64/apps/grass/grass-7.0.5", home = tempdir(),
> > > >  override = TRUE)
> > > >  data(meuse)
> > > >  coordinates(meuse) <- ~x+y
> > > >  writeVECT(meuse, "meuse")
> > > >  execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 
> > > >  double

> > > >  precision")
> > >  This appears to run correctly on Windows 7 with GRASS 7.0.5 Windows
> > >  stand-alone. Do you need to 

Re: [R-sig-Geo] rgrass7 : Error in parseGRASS

2016-11-23 Thread Chris Reudenbach

Michael,


using GRASS from R is a bit cumbersome. GRASS (and this is for sure one 
of it strength ) is extremely rigid with projections etc. Due to this 
you should setup a correct structure. I use a setup that is similar to 
the below snippet.


cheers Chris


library(rgrass7)
library(sp)

# setup a temp workingdir
working.dir<- "~/tmp/"

# get meuse data

data(meuse)
data(meuse.grid)

# georeference the meuse grid data
coordinates(meuse.grid) <- ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
gridded(meuse.grid) <- TRUE

# get a first cellsize/pixel resolution for GRASS
resolution <- meuse.grid@grid@cellsize[1]

# georeference the meuse data
coordinates(meuse) <- ~x+y
proj4string(meuse) <- as.character(CRS("+init=epsg:28992"))

# get projection, proj4 string and extent for GRASS
projection<-(strsplit(meuse@proj4string@projargs,split = " "))
proj4<- paste(projection[[1]][2:length(unlist(projection))], collapse = ' ')
xmax<-meuse@bbox[3]
xmin<-meuse@bbox[1]
ymax<-meuse@bbox[4]
ymin<-meuse@@bbox[2]

# create and set working directory
if (!file.exists(file.path(working.dir,"run"))){
  dir.create(file.path(working.dir,"run"),recursive = TRUE)
}
setwd(file.path( working.dir,"run"))


 set up GRASS RUNTIME environment

# define the GRASS executable path
if(Sys.info()["sysname"] == "Windows"){
  grass.gis.base<-'C:/OSGeo4W64/apps/grass/grass-7.0.5'
}else {
  grass.gis.base<-'/usr/lib/grass72'
}

# set path for optional GRASS addons
Sys.setenv(GRASS_ADDON_PATH="~/.grass7/addons")

# create the TEMPORARY GRASS location
rgrass7::initGRASS(gisBase=grass.gis.base,
   home=tempdir(),
   mapset='PERMANENT',
   override=TRUE
)

# assign GRASS projection according to data set
rgrass7::execGRASS('g.proj',
   flags=c('c','quiet'),
   proj4=proj4
)

# assign GRASS extent and resolution
rgrass7::execGRASS('g.region',
   flags=c('quiet'),
   n=as.character(ymax),
   s=as.character(ymin),
   e=as.character(xmax),
   w=as.character(xmin),
   res=as.character(resolution)
)


#   now do GRASS STUFF

writeVECT(meuse,"meuse")
execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 double
precision")


On 23.11.2016 07:38, Michael DELORME wrote:

Thanks for your insight.
I'll use a stand-alone GRASS install.

Cordially

Le 22/11/2016 15:13, Roger Bivand a écrit :

On Tue, 22 Nov 2016, Michael DELORME wrote:


Thanks for investigating.

I'll try a standalone version of GRASS if needed.

Here is the traceback


traceback()

4: stop(paste(cmd, "not parsed"))
3: parseGRASS(cmd, legacyExec = legacyExec)
2: doGRASS(cmd, flags = flags, ..., parameters = parameters, echoCmd =
echoCmd, legacyExec = legacyExec)
1: execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 double
precision")

My C:\OSGeo4W64\apps\grass\grass-7.0.5\bin lists *.exe or *.bat and
indeed the *.bat command are those which don't work (r.mask.bat,
v.db.addcolumn.bat,...). The batch calls Python scripts ; for example
v.db.addcolumn.bat is :

@"%GRASS_PYTHON%" "%GISBASE%/scripts/v.db.addcolumn.py" %*

So I guess something must be wrong in my Python installation... (it
works however from within a GRASS shell)

No, the logic for Windows and initGRASS has only ever been tested with
stand-alone GRASS. I would guess that in OSGeo4W the path variables
are wrong. I cannot even find out how to run initGRASS in R but
outside the OSGeo4W shell (doesn't find iconv.dll) or inside the
OSGeo4W shell (cannot create a temporary GRASS location). Unless you
or others really need OSGeo4W, I suggest asking others to implement
that - very bulky and fragile.

Roger


Cordially


Le 22/11/2016 14:30, Roger Bivand a écrit :

On Tue, 22 Nov 2016, Michael DELORME wrote:


Thanks for replying

Here is a reproducible example (change your GRASS directory of
course) :

library(rgrass7)
library(sp)
initGRASS("C:/OSGeo4W64/apps/grass/grass-7.0.5", home = tempdir(),
override = TRUE)
data(meuse)
coordinates(meuse) <- ~x+y
writeVECT(meuse, "meuse")
execGRASS("v.db.addcolumn", map = "meuse", columns = "area_m2 double
precision")

This appears to run correctly on Windows 7 with GRASS 7.0.5 Windows
stand-alone. Do you need to use OSGeo4W - the internal logic used for
finding out where the different versions keep their executables and
batch files varies?

What does traceback() say after the error? The error comes from
parseGRASS("v.db.addcolumn") - the XML error comes from not being able
to read the output of v.db.addcolumn<.ext> --interface-description
where <.ext> may be .exe or .bat I think.

Roger


I get :
Error : XML content does not seem to be XML: 'The specified path
was not
found.'
In addition: Warning message:
running command 'v.db.addcolumn.bat --interface-description' had
status 1
Error in parseGRASS(cmd, legacyExec = legacyExec) :
  v.db.addcolumn not parsed

Other execGRASS