Re: [R-sig-Geo] Problem with initGRASS

2018-11-23 Thread Chris Reudenbach 
Jaime,

if the location of your GRASS installation is not known you may also 
give the link2GI package a try.

try link2GI::linkGRASS7() for a "brute force" searching and linking. 
This should setup your enviroment for using the newest GRASS version.

cheers Chris

On 23/11/2018 09:44, Roger Bivand wrote:
> On Fri, 23 Nov 2018, Jaime Burbano Girón wrote:
>
>> Hi everyone,
>>
>> I'm trying to use GRASS in R, but I can't get to initialize it. I'm 
>> using
>> Ubuntu 18.04 and GRASS 7.4.2. This is the system information:
>>
>> Versión de GRASS: 7.4.2
>>
>> Revisión SVN de GRASS: exported
>>
>> Fecha de compilación: 2018-10-28
>>
>> Construir plataforma: x86_64-pc-linux-gnu
>>
>> GDAL: 2.2.3
>>
>> PROJ.4: 4.9.3
>>
>> GEOS: 3.6.2
>>
>> SQLite: 3.22.0
>>
>> Python: 2.7.15rc1
>>
>> wxPython: 3.0.2.0
>>
>> Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic
>>
>>
>> I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
>> And this is the error:
>
> No, ?initGRASS says:
>
> gisBase: The directory path to GRASS binaries and libraries
>
> and grass74/bin only contains (for me):
>
> $ ls bin
> grass74
>
> which is a Python script, ASCII text executable. Trying
>
> grass74/grass-7.4.2
>
> (no idea where your installer puts this):
>
> $ ls grass-7.4.2
> AUTHORS   contributors_extra.csv  fonts REQUIREMENTS.html
> bin   COPYING GPL.TXT  scripts
> CHANGES   demolocation    gui  share
> CITING    docs    include  tools
> config.status driver  INSTALL  translators.csv
> contributors.csv  etc lib
>
> which contains the binaries in bin and the libraries in lib. The 
> required string is what is set in the shell variable GISBASE when 
> running GRASS interactively:
>
> echo $GISBASE
>
> hence the name of the argument. initGRASS() cannot discover this as 
> GRASS is not running, and has not set its environment variables. I 
> have added a test looking for a directory called bin in the directory 
> given as the gisBase argument, check on 
> https://r-forge.r-project.org/R/?group_id=2020 to see when rgrass7 
> 0.1-13 (Rev: 68) has build status current, try that with your current 
> setting, the correct setting and report back.
>
> Hope this clarifies,
>
> Roger
>
>
>>
>> sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
>> not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
>> g.version: not foundError in system(paste("g.version", get("addEXE",
>> envir = .GRASS_CACHE),  :
>>  error in running commandAdemás: Warning messages:1: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command2: In system(paste(paste("g.gisenv",
>> get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command4: In system(paste(paste("g.gisenv",
>> get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
>> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
>> error in running command
>>
>>
>> I verified that g.gisenv and g.version are located in 
>> /usr/lib/grass74/bin.
>> Any suggestion?
>>
>> Thanks in advance.
>>
>> Best,
>>
>> Jaime.
>>
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, 
GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: 
++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, 
moc.environmentalinformatics-marburg.de


[[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] Is the raster package still receiving updates?

2017-10-21 Thread Chris Reudenbach 
Many of us rely on the raster package as a crucial basis for their work. 
Big thank to Robert.  Nevertheless to contribute needs to focus needs 
and ressources.  Maybe an information exchange/discussion with Robert 
about options such as if/how the interested community can support would 
be more productive.


cheers Chris

Am 21.10.2017 um 06:52 schrieb Roger Bivand:

The CRAN raster package has not been abandoned, if it had, it would have been 
reclassified as orphaned. Its test results have notes but no warnings or 
errors. It would be orphaned if the maintainer did not respond to CRAN requests 
to resolve test errors. Consequently, you are suggesting a fork, and should 
rename any package. Only the maintainer may update existing CRAN packages. I 
would not think that forking an existing CRAN package is the most productive 
way of contributing to the community. The package has many reverse 
dependencies, so setting up a test framework to find backwards 
incompatibilities would be crucial.

Of course, all contributions are welcome.

Roger

Roger Bivand
Norwegian School of Economics
Bergen, Norway



Fra: Michael Sumner
Sendt: lørdag 21. oktober, 06.28
Emne: Re: [R-sig-Geo] Is the raster package still receiving updates?
Til: Hodgess, Erin
Kopi: R-sig-geo Mailing List


Great! I'm keen, but haven't done the work needed to flesh this out properly: https://github.com/mdsumner/raster-rforge/issues I have not checked if there's been any commits to r-forge since I cloned this. Cheers, Mike On Sat, 21 Oct 2017, 12:58 Hodgess, Erin, wrote: > I could take a swat at it, if you wish. Please do send me the list of 
> fixes. I could start in about a week, if that would work. > Sincerely, > Erin > > > Erin M. Hodgess > Associate Professor > Department of Mathematics and Statistics > University of Houston - Downtown > mailto: hodge...@uhd.edu > >  > From: R-sig-Geo on behalf of 
Michael > Sumner > Sent: Friday, October 20, 2017 8:05 PM > To: Thiago V. dos Santos > Cc: R-sig-geo Mailing List > Subject: Re: [R-sig-Geo] Is the raster package still receiving updates? > > No, it's effectively abandoned. No concrete plans known. > > I'm interested to keep it going and have kept a list of fixes 
needed, but > not sure when or if I'll get to it. > > I'd support any efforts in this, I'll rely on raster for many years still. > > Cheers, Mike > > On Sat, 21 Oct 2017, 02:05 Thiago V. dos Santos via R-sig-Geo, < > r-sig-geo@r-project.org> wrote: > > > Dear list, > > > > I realized that the 
latest version of the raster package was released on > > CRAN over a year ago. > > > > I also noticed that Robert's participation in this list has become rather > > scarce. > > > > I was just wondering whether the raster package is still receiving > updates? > > > > Greetings, > > -- 
Thiago V. dos Santos > > > > Postdoctoral Research Fellow > > Department of Climate and Space Science and Engineering > > University of Michigan > > > > ___ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > 
https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > -- > Dr. Michael Sumner > Software and Database Engineer > Australian Antarctic Division > 203 Channel Highway > Kingston Tasmania 7050 Australia > > [[alternative HTML version deleted]] > > 
___ > R-sig-Geo mailing list > 
R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Dr. 
Michael Sumner Software and Database Engineer Australian Antarctic Division 203 Channel 
Highway Kingston Tasmania 7050 Australia [[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



--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, 
GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: 
++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, 
moc.environmentalinformatics-marburg.de

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

Re: [R-sig-Geo] Optimizing spatial query in R

2017-02-20 Thread Chris Reudenbach 

Ervan,

Just on the run, especially if you deal with real data and some 100K 
-1000K of vectors I think the fastest way to do so is to engage GDAL 
GRASS7/SAGA and their ability to deal highly efficient with this kind 
of topological and geometrical queries.


A very good pure R alternative is to use the sf package that is ways 
faster and more satble than sp. It also provides this type of GIS 
functionality like st_overlaps() etc.


cheers Chris



On 20.02.2017 09:45, Ervan Rutishauser wrote:

Dear All,

I am desperately trying to fasten my algorithm to estimate the fraction of
tree crown that overlap a given 10x10 subplot in a forest plot. I have
combined a set of spatial functions (gDistance, extract)  & objects
(SpatialGrid, SpatialPolygons) in a way that is probably not the most
efficient, as it takes dozen of hours to run for a 5 subplots (50 ha
forest plot).

My detailed problem and a reproducible example are posted on Stackoverflow
 and
append below, if you wanna have a look. Apart of Amazon Web Server, is
anyone aware of a sever where to execute (and save results back) R codes
online?
Thanks for any help.
Best regards,


# A. Define objects
 require(sp)
 require(raster)
 require(rgdal)
 require(rgeos)
 require(dismo)
 radius=25   # max search radius around 10 x 10 m cells
 res <- vector() # where to store results

 # Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)
 set.seed(0)
 survey <- 
data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))
 coordinates(survey) <- ~x+y

 # Define 10 x 10 subplots
 grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10)))
 survey$subplot <- over(survey,grid10)
# B. Now find fraction of tree crown overlapping each subplot
 for (i in 1:100) {
 # Extract centroïd of each the ith cell
 centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]
 corner <-
data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))


 # Find trees in a max radius (define above)
 tem <- 
survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]


 # Define tree crown based on tree diameter
 tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

 # Compute the distance from each tree to cell's borders
 pDist <- vector()
 for (k in 1:nrow(tem))  {
 pDist[k] <-
gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1
 }

 # Keeps only trees whose crown is lower than the above
distance (=overlap)
 overlap.trees <- tem[which(pDist<=tem$crownr),]
 overlap.trees$crowna <-overlap.trees$crownr^2*pi  # compute crown area

 # Creat polygons from overlapping crowns
 c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr,
lonlat=F, dissolve=F)
 crown <- polygons(c1)
 Crown <-
SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

 # Create a fine grid points to retrieve the fraction of
overlapping crowns
 max.dist <- ceiling(sqrt(which.max((centro$x -
overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance
to narrow search

 finegrid <-
as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))
 coordinates(finegrid) <- ~ x+y
 A <- extract(Crown,finegrid)
 Crown@data$ID <- seq(1,length(crown),1)
 B <- as.data.frame(table(A$poly.ID))
 B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)
 B$overlap <- B$Freq/B$crown.area
 B$overlap[B$overlap>1] <- 1
 res[i] <- sum(B$overlap)
 }
# C. Check the result
 res  # sum of crown fraction overlapping each cell (works fine)



--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of 
Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 
Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: 
gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de


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

Re: [R-sig-Geo] agent-based model with GIS

2017-02-03 Thread Chris Reudenbach 

Antony,

with focus on prototyping and medium to complex modeling I strongly 
suggest Uri Wilensky's  Netlogo (https://ccl.northwestern.edu/netlogo/) 
which can interfaced from R using RNetLogo 
https://cran.r-project.org/web/packages/RNetLogo/index.html that 
provides an excellent wrapper .


cheers Chris

On 03.02.2017 16:43, Antony wrote:

Dear all,

I am interested in building an agent-based model with geospatial information in R. 
My idea is to create agents that walk down streets, drive cars and generally act 
like ordinary people going to work and so on. I was wondering if this is something 
that anyone on this group has tried to do using R or knows about? Specifically, I 
am interested in using this to understand how crowds form. So far, I have been 
building simulations like this in Netlogo, but am interested to move into using R 
and would very much appreciate people's experience / hints & tips.

Many thanks & best wishes,

Antony

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



--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, 
GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: 
++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, 
moc.environmentalinformatics-marburg.de

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

Re: [R-sig-Geo] Distance calculating between points according to a river network

2017-01-25 Thread Chris Reudenbach 

Tristan,

if I understand you right you just want to calculate the distances. So 
no hydro stuff like accumulation flow directions etc.? As well as no 
diffusion modeling of the pollutants.


Most of the traditional vector based GIS stuff is packed into the rgeos 
package. So you might calculated the distances along a graph like in the 
following paste and copy snippet.


cheers Chris



library(rgeos)
library(sp)
# just for visualisation
r <- raster(nrows=6, ncols=7, xmn=0, xmx=7, ymn=0, ymx=6, crs="+proj=utm 
+units=m")

r[] <- c(1, 0, 1, 0, 1, 0, 1,
 0, 1, 0, 1, 0, 1, 0,
 1, 0, 1, 0, 1, 0, 1,
 0, 1, 0, 1, 0, 1, 0,
 1, 0, 1, 0, 1, 0, 1,
 0, 1, 0, 1, 0, 1, 0)
# river line
line 
<-SpatialLines(list(Lines(Line(cbind(c(5.5,1.5,4.5),c(1.5,5.5,5.5))), 
ID="line")))

# points
p1 <-as.data.frame(cbind(4.5,2.5))
p2 <-as.data.frame(cbind(2.5,4.5))
p3 <-as.data.frame(cbind(2.1,4.5))
p4 <-as.data.frame(cbind(3.5,5.5))
p5 <-as.data.frame(cbind(3.5,5.7))
sp::coordinates(p1) <- ~V1+V2
sp::coordinates(p2) <- ~V1+V2
sp::coordinates(p3) <- ~V1+V2
sp::coordinates(p4) <- ~V1+V2
sp::coordinates(p5) <- ~V1+V2
projection(line)<-projection(point)

plot(r)
plot(line,add=TRUE)
plot(p1,add=TRUE)
plot(p2,add=TRUE)
plot(p3,add=TRUE,col="red")
plot(p4,add=TRUE)
plot(p5,add=TRUE,col="red")

# calculate distance along the line
rgeos::gProject(spgeom = line, sppoint = p1)
rgeos::gProject(spgeom = line, sppoint = p2)
rgeos::gProject(spgeom = line, sppoint = p4)

# note if the point is not near the line it is caluclated orthogonal
rgeos::gProject(spgeom = line, sppoint = p3)
rgeos::gProject(spgeom = line, sppoint = p5)



On 25.01.2017 17:31, Romine, Jason wrote:

The gdistance package in R is what you are looking for:

https://cran.r-project.org/web/packages/gdistance/gdistance.pdf

https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf

Best,
Jason

On Wed, Jan 25, 2017 at 5:22 AM, Tristan Bourgeois <
tristan.bourge...@gmail.com> wrote:


Hi Mirza,

Thanks for your quick answer.

Actually I don't want to calculate distance between a point and the next
downstream one.
For example, if there are 4 sewage plant on a river and 7 measurment
stations,I want the distance between each sewage plant and each station.
(water pollution dispersion calculation).
So the flow direction is not that important in my case. The difficulty is
to get the distance between a point and another one following the river
polyline.

In your case, you can add a raster layer with elevation values ans extract
the elevation  for each point of your layer.  Then, for each polyline it is
possible to define a flow direction.
Hope it will help you.




2017-01-25 11:57 GMT+01:00 Mirza Cengic :


Hi Tristan,

I recently had to do something similar (calculating distance to the next
downstream point), and I was not able to find a solution with R that

would

deal with the direction of the flow. Perhaps someone who has deeper
understanding of packages such as iGraph etc. might be able to find a
solution, but I wasn't.
I solved the problem with ArcMap, using HydroTools extension, where you
need to define a geometric network that has a direction, and then the
calculation is pretty straightforward with HydroTools (let me know if you
need more info).
I would prefer to have a scripted solution in an open-source software,

but

this worked for me in any case. If anyone has a better way for solving
this, I would love to see how!

Cheers,
Mirza.




--
Tristan Bourgeois

 [[alternative HTML version deleted]]

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






--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, 
GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: 
++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, 
moc.environmentalinformatics-marburg.de

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

Re: [R-sig-Geo] shortestPath in a loop

2017-01-25 Thread Chris Reudenbach 
Marta,

your problem seems to me more conceptual. The shortestPath 
implementation of a cost analysis is not exactly want you want. You 
would choose it if you have to find an unknown path on a friction 
surface. You would not choose it if you just want to link positions. The 
shortestPath of your implementation can not avoid to touch the island 
because you force it to go there. It is a bit more "expensive" in your 
case 2 units but not forbidden.

Nevertheless there are also some technical issues to address: Even if if 
there is a correction for unprojected data sets it is better to project 
them because it is more stable to calculate cost on an equal area 
projection. Than the spatial resolution of the cost raster should be 
sufficient to meet the distance of the points.

You will find below the example code. Note that I assume that the data 
is an directory as defined by "path_data".

The "path will be plotted consecutively in the plotting window.

hope this clarifies

cheers Chris

#question 4
library(sp)
library(maptools)
library(rgeos)
library(mapview)
library(robubu)
library(raster)

### define data folder
path_data<-"~/proj/boats/data/"

### boat points the locations of the boats are sort by date and time!
boat <- read.table(paste0(path_data,"boat2905.csv"), header=TRUE, 
sep=",", na.strings="NA", dec=".", strip.white=TRUE)
pos<-as.data.frame(cbind(boat$Lat1,boat$Long1))
sp::coordinates(pos) <- ~V2+V1
sp::proj4string(pos) <-CRS("+proj=longlat +datum=WGS84 +no_defs")

### project it to somewhat equal area
boatSP <- sp::spTransform(pos, CRS("+proj=laea +lat_0=30 +lon_0=-30 
+x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))

# cost raster sea and island
azoTS1<- raster::raster(paste0(path_data,"azoTS1.tif"))

### project it to somewhat equal area
azoTS1<- raster::projectRaster(azoTS1,  crs="+proj=laea +lat_0=30 
+lon_0=-30 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")

###scale  real cost file -> makes it MUCH more expansive to sail over 
land surfaces
azoTS1[azoTS1 ==2]<-100

### increase the spatial resolution
costraw<-disaggregate(azoTS1,fact=c(10, 10))

# prepare it for analysis
transition<- gdistance::transition(1/costraw, mean, directions=8)
trCostS16 <- gdistance::geoCorrection(transition, type="c")

### run first analysis -> will link all boat positions
### note it will NOT avoid sailing over land because it is forced to do 
so by setting the last position on the (rasterized) island
raster::plot(azoTS1)
pathwaylist<- lapply(seq(2:length(boatSP)), function(i){
   costpath<- try(shortestPath(trCostS16, 
boatSP@coords[i,],boatSP@coords[i-1,], output="SpatialLines"),silent = TRUE)
   lines(costpath)
})

### now we import the species point
### just for showing that the cost analysis will work
spcs <- read.table(paste0(path_data,"/sp2pontosDT.csv"), header=TRUE, 
sep=",", na.strings="NA", dec=".", strip.white=TRUE)#Bom test
sp::coordinates(spcs) <- ~Long1+Lat
sp::proj4string(spcs) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
spcsSP <- sp::spTransform(spcs, CRS("+proj=laea +lat_0=30 +lon_0=-30 
+x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))

### run second analysis -> will link all the two species positions
raster::plot(azoTS1)
pathwaylist<- lapply(seq(2:length(spcsSP)), function(i){
   costpath<- try(shortestPath(trCostS16, 
spcsSP@coords[i,],spcsSP@coords[i-1,], output="SpatialLines"),silent = TRUE)
   lines(costpath)
})



On 24.01.2017 21:06, marta azores wrote:
> Hi guys,
>
> I'm still trying to solve my deal. How I can join the survey 
> boat points using the shortestpath function, to avoid the land. It's 
> important because  I need an output "spatialLines"(boat surveys) to 
> intersect with points( species points).
> I add several files : script and other layers of information.
>
> Thanks in advance,
>
> Marta
>
> 2017-01-16 9:45 GMT-01:00 marta azores  >:
>
> Hi Jerome,
>
> There are several R packages to join points. However, I have two
> relevant issues with them.
>  -The first,  I need to get the "Lines" from these tracks in
> "spatial Lines". For example, in the adehabitatLT package you can
> get easily the trajectories of the animal but with the class "
> ltraj". I don't know
>how transform this object class into spatial lines.
> -The second, when I made the tracks I need to avoid the coast, and
> I don't know how do that in adehabitatLT.
>
> Any idea?
>
> Thanks for your suggestions,
>
> Marta
>
> 2017-01-14 15:21 GMT-01:00 Roland Harhoff
>  >:
>
> Hi Marta,
>
> and did you check the trajectories package ...?
>
> https://cran.r-project.org/web/packages/trajectories/index.html
> 
>
> Cheers!
> Roland
>
>
> J�rome Mathieu schrieb am 2017-01-14:
>

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

2016-11-24 Thread Chris Reudenbach 
@Michael fine!

@Roger

this is  pretty dirty "hack" I did not cross check if I got all 
variables dealing with GRASS. I would like to do it more systematically. 
Nevertheless as far as I understand the osgeo4w concept, we have only 
two changing names. 1) the installation root directory   2) the GRASS 
version name, that is used as a part of the internal  path structure.  I 
would like to check this and if so it might be easy to integrate this in 
initGRASS.

cheers Chris

Am 24.11.2016 um 10:48 schrieb Roger Bivand:
> On Thu, 24 Nov 2016, Michael DELORME wrote:
>
>> That's great Chris, it works fine as is !
>> Thank you, I don't have to keep another GRASS install aside then.
>
> Thanks for confirming this. I'll try to set up a test rig to allow 
> initGRASS to set up the environment variables when using GRASS in 
> OSGeo4W for a later release.
>
> Roger
>
>>
>> I'm sure it will be helpful to other people also.
>> Thanks again
>>
>> Le 24/11/2016 09:17, Chris Reudenbach a �crit :
>>> Michael,
>>>
>>> Using a fresh OSGEO4W64 standard quick desktop installation (i.e.
>>> installation of the default desktop GIS software packages with the
>>> installer and using C:\OSGeo4W64 as path)  I suggest the below 
>>> solution.
>>>
>>> You have to set all necessary eniromental and system variables
>>> manually. I am not quit sure if I got them but actually it seems to
>>> work for your question. As a conclusion I think your Python and all of
>>> the rest is installed correct.
>>>
>>>
>>> cheers Chris
>>>
>>>
>>> ### librarys, data etc
>>> -
>>> 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) <- 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"))
>>>
>>>
>>> ### SETUP OSGEO4W enviroment settings manually
>>> -
>>> # setup the OSGEO4W environ manually
>>> # assuming a osgeow4w default "deskop fastinstall
>>> # using the default installation directory "C:\OSGeo4W64"
>>>
>>> # set OSGE4W base directory
>>> osgeo4w.root<-"C:\\OSGEO4~1"
>>> Sys.setenv(OSGEO4W_ROOT=osgeo4w.root)
>>> # define GISBASE
>>> grass.gis.base<-paste0(osgeo4w.root,"\\apps\\grass\\grass-7.0.5")
>>> Sys.setenv(GISBASE=grass.gis.base)
>>>
>>> Sys.setenv(GRASS_PYTHON=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\bin\\python.exe"))
>>>  
>>>
>>>
>>> Sys.setenv(PYTHONHOME=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\apps\\Python27"))
>>>  
>>>
>>>
>>> Sys.setenv(PYTHONPATH=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\apps\\grass\\grass-7.0.5\\etc\\python"))
>>>  
>>>
>>>
>>> Sys.setenv(GRASS_PROJSHARE=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\proj"))
>>>  
>>>
>>>
>>> Sys.setenv(PROJ_LIB=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\proj"))
>>> Sys.setenv(GDAL_DATA=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\gdal")) 
>>>
>>> Sys.setenv(GEOTIFF_CSV=paste0(Sys.getenv("OSGEO4W_ROOT"),"

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

2016-11-24 Thread Chris Reudenbach 

Michael,

Using a fresh OSGEO4W64 standard quick desktop installation (i.e. 
installation of the default desktop GIS software packages with the 
installer and using C:\OSGeo4W64 as path)  I suggest the below solution.


You have to set all necessary eniromental and system variables manually. 
I am not quit sure if I got them but actually it seems to work for your 
question. As a conclusion I think your Python and all of the rest is 
installed correct.



cheers Chris


### librarys, data etc 
-

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


### SETUP OSGEO4W enviroment settings manually 
-

# setup the OSGEO4W environ manually
# assuming a osgeow4w default "deskop fastinstall
# using the default installation directory "C:\OSGeo4W64"

# set OSGE4W base directory
osgeo4w.root<-"C:\\OSGEO4~1"
Sys.setenv(OSGEO4W_ROOT=osgeo4w.root)
# define GISBASE
grass.gis.base<-paste0(osgeo4w.root,"\\apps\\grass\\grass-7.0.5")
Sys.setenv(GISBASE=grass.gis.base)

Sys.setenv(GRASS_PYTHON=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\bin\\python.exe"))
Sys.setenv(PYTHONHOME=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\apps\\Python27"))
Sys.setenv(PYTHONPATH=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\apps\\grass\\grass-7.0.5\\etc\\python"))
Sys.setenv(GRASS_PROJSHARE=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\proj"))
Sys.setenv(PROJ_LIB=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\proj"))
Sys.setenv(GDAL_DATA=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\gdal"))
Sys.setenv(GEOTIFF_CSV=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\epsg_csv"))
Sys.setenv(FONTCONFIG_FILE=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\etc\\fonts.conf"))

# call all OSGEO4W settings
system("C:/OSGeo4W64/bin/o-help.bat")


# create PATH variable
Sys.setenv(PATH=paste0(grass.gis.base,";",
"C:\\OSGEO4~1\\apps\\Python27\\lib\\site-packages\\numpy\\core",";",
"C:\\OSGeo4W64\\apps\\grass\\grass-7.0.5\\bin",";",
"C:\\OSGeo4W64\\apps\\grass\\grass-7.0.5\\lib",";",
"C:\\OSGeo4W64\\apps\\grass\\grass-7.0.5\\etc",";",
"C:\\OSGeo4W64\\apps\\grass\\grass-7.0.5\\etc\\python",";",
"C:\\OSGeo4W64\\apps\\Python27\\Scripts",";",
   "C:\\OSGeo4W64\\bin",";",
   "c:\\OSGeo4W64\\apps",";",
   "C:\\OSGEO4~1\\apps\\saga",";",
   paste0(Sys.getenv("WINDIR"),"/WBem"),";",
   Sys.getenv("PATH")))


 start with GRASS setup 


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 
-


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

   precision")
#  do other command line stuff 
-


## call gdal
system("gdal_merge")

# call saga cli
system("saga_cmd")


Am 23.11.2016 um 07:38 schrieb Michael DELORME:

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

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

Re: [R-sig-Geo] best practice for reading large shapefiles?

2016-04-26 Thread Chris Reudenbach 

Vinh

Even if it might be in this list OT, IMHO R is not the best tool for 
dealing with this amount of vector data. Actually I agree completely 
with Roger's remarks and corresponding to the "competent platform" you 
also may think about using software for big data...


As Roger already has clarified: The recommendation what might be best 
depends highly  on your questions and issues or on the type of analysis 
you need to run and cannot be answered straightforward.


I think Edzer can clarify up to which size sp object are still "usable", 
following my experience  i would guess something like 500K polygons 1M 
lines and up to 5M points but it is highly dependent on the number of 
attributes. So you are far beyond this.



If you want to deal with this amount of spatial vector data using R, it 
is highly reasonable to have a look at one of the mature GIS packages 
like GRASS or QGIS. You can use them via their APIs.
Nevertheless you easily can put it in postgres/postgis and perform all 
operations/analysis using the spatial capabilities and build in 
functions of postgis if you are an experienced PostGis user.


cheers
Chris


Am 26.04.2016 um 22:33 schrieb Vinh Nguyen:

On Tue, Apr 26, 2016 at 1:12 PM, Roger Bivand  wrote:

On Tue, 26 Apr 2016, Vinh Nguyen wrote:


Would loading the shapefile into postgresql first and then use readOGR
to read from postgres be a recommended approach?  That is, would the
bottleneck still occur?  Thank you.


Most likely, as both use the respective OGR drivers. With data this size,
you'll need a competent platform (probably Linux, say 128GB RAM) as
everything is in memory. I find it hard to grasp what the point of doing
this might be - visualization won't work as none of the considerable detail
certainly in these files will be visible. Can you put the lot into an SQLite
file and access the attributes as SQL queries? I don't see the analysis or
statistics here.


- I can't tell from your response whether you are recommending PostGIS
is a recommended approach or not.  Could you clarify?

- I am working on a Windows server with 64gb ram, so not too weak,
especially for some files that are a few gb in size.  Again, not sure
if the job just halted or it's still running, but just rather slow.
I've killed it for now as the memory usage still has not grown after a
few hours.

- Yes, the shapes are quite granular and many in quantity.  The use
case was not to visualize them all at once.  Wanted a master file so
that when I get a data set of interest, I could intersect the two and
then subset the areas of interest (eg, within a state or county).
Then visualize/analyze from there.  The master shapefile was meant to
make it easy (reading in one file) as opposed to deciding which
shapefile to read in depending on the project.

- I just looked back at the 30 PLSS zip files, and they provide shapes
for 3 levels of granularity.  I went with the smallest.  I just
realized that the mid-size one would be sufficient for now, which
results in dbf=138mb and shp=501mb.  Attempting to read this in now (~
30 minutes), which I assume will read in fine after some time.  Will
respond to this thread if this is not the case.

Thanks for responding Roger.

-- Vinh

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




--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of 
Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 
Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: 
gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de


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

Re: [R-sig-Geo] Find a circle center with spatial points

2016-03-19 Thread Chris Reudenbach 
Because it seems to be an arc and not a circle issue that you can solve 
the problem by
picking arbitrary two points of your assumed "arc" then construct 
(calculate)  the perpendicular bisector of
the line between them and do so for another arbitrary two points of the 
assumed "arc".


The intersection of the perpendicular lines is the assumed center of the 
arc.


If you iterate over all points this should be a pretty good estimation 
of the real center.


cheers Chris

Am 18.03.2016 um 18:36 schrieb Barry Rowlingson:

On Fri, Mar 18, 2016 at 2:43 PM, Alex Mandel  wrote:

library(rgeos)
gCentroid

http://www.rdocumentation.org/packages/rgeos/functions/topo-unary-gCentroid

Assuming its a circle that would be the center.

Only if you have points uniformly (or uniform-randomly) distributed
round the full extent of the circle. From Adrien's plot it looks like
he's got an arc there.

  It seems more like a three-parameter optimisation problem. Find x, y,
and r that define the circle that minimises the sum of squared
distances from data points to the circle.

  I'm not sure how you'd choose a good initial x,y,r for your optimiser
since I suspect the surface you're optimising over is not unimodal...
You could try taking lots of random samples of three points from your
data and computing the unique circle that fits those points, then
using the mean (or possibly median, there's a fair chance of massive
outliers) value as the initial values.

A quick googling has actually found this little paper on the subject:

http://www.spaceroots.org/documents/circle/circle-fitting.pdf

So I'll shut up now.

Barry

___
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] Find a circle center with spatial points

2016-03-18 Thread Chris Reudenbach 

Oops :-[ , Barry thanks for clarification , I got it now.

I don't understand the math either but as a kind of "reparation"
for not reading well enough I found the authors java code.
https://www.spaceroots.org/documents/circle/CircleFitter.java
After download you should for standalone compilation delete line 37 
"package org.spaceroots;"

compile it with:

 javac CircleFitter.java -Xlint:unchecked

then you can run it with:

java CircleFitter input.file

The input data should be formated like:

# input
# x y
0 5
1 4.5
2.5 4
3 3.5
4 2
5 0

the above input yields

initial circle: -001.69467803 -000.69446643 006.23578103
converged after 7 iterations
final circle: -001.34339845 -001.34426151 006.44308386
with the format x,y,radius

which at least make sense with respect to the data.

You easily can  run it from R by system() and grab the output in a file 
using a OS depending pipe

or you may run it using the rJava package which could be more complex.

cheers Chris

if you want you can call it from R with javaR

Am 18.03.2016 um 20:11 schrieb Barry Rowlingson:

On Fri, Mar 18, 2016 at 6:49 PM, Chris Reudenbach
<reudenb...@uni-marburg.de> wrote:

Because it seems to be an arc and not a circle issue that you can solve the
problem by
picking arbitrary two points of your assumed "arc" then construct
(calculate)  the perpendicular bisector of
the line between them and do so for another arbitrary two points of the
assumed "arc".

The intersection of the perpendicular lines is the assumed center of the
arc.

If you iterate over all points this should be a pretty good estimation of
the real center.

  This is the  "sample 3 points and find the fitted circle" idea, you
are likely to get massive "outliers" and if you take the mean
coordinate it could fail horribly. See the paper I linked to for an
example. They use the median to get an initial "robust" estimate of
x,y,R, and then use some specialised optimisation to improve the
estimate - you can't just throw it into "optim"!

  Not sure I understand the maths in it yet though

Barry



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

Re: [R-sig-Geo] resolution of openmap() raster layers

2016-02-26 Thread Chris Reudenbach 

Agus,

Mapview is using leaflet as engine. Due to this you will have the 
control icons on the map because first of all it is designed for 
interactive mapping within RStudio/R.


I think there are two different approaches to save your maps:

If you want to have a dump of the mapviewobject (but including the 
graphical buttons) you'll find a descriptat stackoverflow

http://stackoverflow.com/questions/31336898/how-to-save-leaflet-in-rstudio-map-as-png-or-jpg-file


You may also use  the spplot function from mapview which is designed for 
basic static mapping and to make usable the adavantages of spplot.


Note even if you are dealing with the mapview map object the spplot 
function  uses the Openstreetmap package for retrieving  the background 
maps  (e.g. 
http://www.inside-r.org/packages/cran/OpenStreetMap/docs/openmap). You 
can use the spplot syntax for designing your maps. Up to now this static 
plotting function is still pretty  basic but you may have a try:
spplot(mapView(nica["POP2000"]),colorkey=FALSE,  lwd= 15, alpha.regions 
= 0.9,

map.type="stamen-watercolor" )

I think for using the spplot it is better to install the current stable 
from github:

library(devtools)
install_github("environmentalinformatics-marburg/mapview", ref = "master")


cheers chris


Am 26.02.2016 um 12:27 schrieb Agustin Lobo:

Stunning!
Can I remove the buttons for saving to a bmp file?
What attribution should be used for publishing?
Agus

On Thu, Feb 25, 2016 at 7:42 PM, Chris Reudenbach
<reudenb...@uni-marburg.de> wrote:

Hi,

if you just want to map the data, mapview could be an option that among
others avoid the pixel stretching.

require(mapview)
require(raster)
nica <- getData("GADM", country="NIC", level=0)

mapview(nica)

mapview(nica,zcol = "POP2000", color = "#FFA500", lwd= 5, alpha.regions =
0.4)


cheers Chris




Am 25.02.2016 um 18:49 schrieb Barry Rowlingson:

On Thu, Feb 25, 2016 at 5:11 PM, Agustin Lobo <alobolis...@gmail.com>
wrote:

Is there any way to download the raster layers
of openmap() with an increased resolution?
I find the quality of the labels very low,
or am I doing something wrong? i.e.

require(raster)
require(mapmisc)
nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
plot(nicabg)

   Map tiles from OpenStreetMap and other map tile providers are images
designed to be shown at a fixed resolution. When you plot them in an R
graphics window you could be stretching them so that each pixel in the
original maps to 1.273 pixels on your screen. So some kind of
interpolation or nearest neighbour replacement has to be done, and
this makes text labels look bad. Other line work will look bad too.

   If you try and download more map tiles at a higher resolution then
you'll find the labels are now way too small, because what you've
downloaded are map tiles designed for a higher zoom level on a web
browser. Web map browsers have a fixed set of zoom values that
correspond to the resolution of the map tiles. With an R window, you
are free to choose odd zoom factors that give the ugly behaviour.

   If you can resize your R window exactly right then you might get
something that looks good!

   The alternative is to build a background map yourself from
OpenStreetMap *vector* data and some code and some styling. Or use a
map tile provider that doesn't have text labels and add them to
selected places with R graphics commands. Lines and polygons will
still be stretched and a bit "jaggy" but our eyes don't notice this as
much as badly scaled text.

Barry

___
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 mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] resolution of openmap() raster layers

2016-02-25 Thread Chris Reudenbach 

Hi,

if you just want to map the data, mapview could be an option that among 
others avoid the pixel stretching.


require(mapview)
require(raster)
nica <- getData("GADM", country="NIC", level=0)

mapview(nica)

mapview(nica,zcol = "POP2000", color = "#FFA500", lwd= 5, alpha.regions 
= 0.4)



cheers Chris



Am 25.02.2016 um 18:49 schrieb Barry Rowlingson:

On Thu, Feb 25, 2016 at 5:11 PM, Agustin Lobo  wrote:

Is there any way to download the raster layers
of openmap() with an increased resolution?
I find the quality of the labels very low,
or am I doing something wrong? i.e.

require(raster)
require(mapmisc)
nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
plot(nicabg)

  Map tiles from OpenStreetMap and other map tile providers are images
designed to be shown at a fixed resolution. When you plot them in an R
graphics window you could be stretching them so that each pixel in the
original maps to 1.273 pixels on your screen. So some kind of
interpolation or nearest neighbour replacement has to be done, and
this makes text labels look bad. Other line work will look bad too.

  If you try and download more map tiles at a higher resolution then
you'll find the labels are now way too small, because what you've
downloaded are map tiles designed for a higher zoom level on a web
browser. Web map browsers have a fixed set of zoom values that
correspond to the resolution of the map tiles. With an R window, you
are free to choose odd zoom factors that give the ugly behaviour.

  If you can resize your R window exactly right then you might get
something that looks good!

  The alternative is to build a background map yourself from
OpenStreetMap *vector* data and some code and some styling. Or use a
map tile provider that doesn't have text labels and add them to
selected places with R graphics commands. Lines and polygons will
still be stretched and a bit "jaggy" but our eyes don't notice this as
much as badly scaled text.

Barry

___
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] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

2016-02-24 Thread Chris Reudenbach 

Agus,

sorry for the addon but I think I have to provide a correction of the 
corner coordinates (e.g. the extent values):


In the example that you have posted below I did calculate the extent 
using the domain center coordinate and the WRF grid resolution in meter 
and the number of rows and cols.


Since Dominik provides a link to the file description of the netcdf file 
I think it is more accurate to reproject the  corner coordinates  as 
given by the netcdf header variables (NC_GLOBAL#corner_lats, 
NC_GLOBAL#corner_lons). Assuming that your variable "dctmx" (which I can 
not identfy it in the nc file)  is of type "Mass" (staggered = M) the 
correct corner coordinates are stored as the first 4 entrys of the dump 
snippet below:


lats<-" 
NC_GLOBAL#corner_lats={12.355667,50.26619,50.26619,12.355667,12.308136,50.18787,50.18787,12.308136,12.210785,50.403816,50.403816,12.210785,12.163345,50.325382,50.325382,12.163345}" 

lons<-" 
NC_GLOBAL#corner_lons={-131.43678,-151.29639,-48.703613,-68.563232,-131.5851,-151.51157,-48.488434,-68.414917,-131.38828,-151.41891,-48.581085,-68.611725,-131.53641,-151.63441,-48.36557,-68.463593}"


after cleaning and converting the strings you may calculate the corner 
coordinates:



library(proj4)
## project mass corner coordinates to lambertian
llMass <- ptransform(cbind(clon[1],clat[1])/180*pi,'+proj=longlat 
+datum=WGS84 +no_defs',proj)
ulMass <- ptransform(cbind(clon[2],clat[2])/180*pi,'+proj=longlat 
+datum=WGS84 +no_defs',proj)
lrMass <- ptransform(cbind(clon[3],clat[3])/180*pi,'+proj=longlat 
+datum=WGS84 +no_defs',proj)
urMass <- ptransform(cbind(clon[4],clat[4])/180*pi,'+proj=longlat 
+datum=WGS84 +no_defs',proj)

wrfLccExtMass<-extent(ulMass[1],lrMass[1],llMass[2],ulMass[2])


According to this the correct extent for mass variables should be:

extent(wrfLccExtMass)
class   : Extent
xmin: -3575343
xmax: 3575342
ymin: -2293058
ymax: 2306330


hope this is correct now

cheers Chris

Am 24.02.2016 um 05:12 schrieb Agus Camacho:

With the help of everybody i finally got this to work, so here is a script
that does the job of reprojecting both, a raster layer obtained from a .nc
and some locations in order to overplot them using plot.raster or mapview.
Im using a combination of the advices of Dominik, Michael and Chris.


require(ncdf4)
require(raster)


setwd("C:/~")


r=raster("geo_em.d01.nc",
   varname="dctmx")# days of ctmax events

# Set extent and projections of rasters for plotting
# chris gave me the orig data fom the nc file because i could not install
gdal
xmin= -3545999
xmax= 3546000
ymin = -2286000
ymax=2286000
pr<- "+proj=lcc +lat_1=25 +lat_2=45 +lat_0=38.08 +lon_0=-100 +x_0=0
+y_0=0"

wrfLccExt<-extent(xmin,xmax,ymin,ymax)

extent(r) <- extent(wrfLccExt)
projection(r) <- pr

# get and prepare  urosaurus locations

x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
x=x[,1:3]
colnames(x)
coordinates(x)=~lon+lat
proj4string(x)<- CRS("+proj=longlat +datum=WGS84")
x=spTransform(x, pr)


# get prepare and plot wrld_simpl
require(maptools)
require(sp)
data(wrld_simpl)
namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
"Canada", "Mexico")), prj)

#plot with raster
plot(r, cex.main=.7,legend=F)
points(x)
plot(namer, add = TRUE)

# plot with mapview (cool!)
m=mapview(r)
u=mapview(x)
m+u


Thanks to all!
Agus

2016-02-23 13:11 GMT-07:00 Agus Camacho :


Thanks for that Dominik,

Giving that projection to either the locations, the raster layer generated
from the .nc file, or both, still did not work. I keep having locations
that should be on land falling far on the sea. Might this be a problem
derived from using raster with a file whose original grid distances are not
constant?

Here is a link with the original file which has the original coordinate
data.

https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0


2016-02-23 12:07 GMT-07:00 Dominik Schneider <
dominik.schnei...@colorado.edu>:


This looks like WRF  data. I just
dealt with this.
The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
+a=637 +b=637 +units=m

+proj=lcc which is usually what wrf is run with.
The tricky part is:
+lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
because every WRF run is different (the WRF Preprocessing System
optimizes the projection for the domain).
and then there is probably no shift so you need(?) +x_0=0 +y_0=0

This gives:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
+ellps=sphere +a=637 +b=637 +units=m +no_defs

But, wrf users like to give out lat and  long so you need to assign it:
+proj=longlat +ellps=sphere +a=637 +b=637 +units=m +no_defs

and then reproject the lat/long to lcc coordinates using this string:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
+ellps=sphere +a=637 +b=637

Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

2016-02-23 Thread Chris Reudenbach 

Agus

sorry I did miss the crucial part:

if you are doing as suggestest you have to define manually the 
Lambertian ymin xmax ymin ymax using the header information of the nc file.


here an example for the unstaggered  data

library(gdalUtils)
library(raster)
library(proj4)
library(mapview)

r<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
  of="GTiff",
  ot="Byte",
  output_Raster=TRUE,
  verbose=TRUE)

finfo <- gdalinfo("geo_em.d01.nc")
##extract parameters
lat_1=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT1", 
finfo))],"=")[[1]][2])
lat_2=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT2", 
finfo))],"=")[[1]][2])
lat0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LAT", 
finfo))],"=")[[1]][2])
lon0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LON", 
finfo))],"=")[[1]][2])
dx=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DX", 
finfo))],"=")[[1]][2])
dy=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DY", 
finfo))],"=")[[1]][2])
y=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG", 
finfo))],"=")[[1]][2])
x=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG", 
finfo))],"=")[[1]][2])

x_0=0
y_0=0

## generate compliant prj4 string
projLcc=paste("+proj=lcc +lat_1=",lat_1," +lat_2=",lat_2,
   " +lat_0=",lat0," +lon_0=",lon0, " +x_0=",x_0," 
+y_0=",y_0,sep="")


## project centre coordinates
tr <- ptransform(cbind(lon0, lat0)/180*pi,'+proj=longlat +datum=WGS84 
+no_defs',proj)


## calculate extent using the Lambertian units (m)
xmin=as.integer(tr[1,1]-(x/2*dx-dx))
xmax=as.integer(tr[1,1]+(x/2*dx-dx))
ymin=as.integer(tr[1,2]-(y/2*dy-dy))
ymax=as.integer(tr[1,2]+(y/2*dy-dy))
wrfLccExt<-extent(xmin,xmax,ymin,ymax)
extent(r) <- extent(wrfLccExt)
projection(r) <- projLcc
mapview(r)


cheers Chris


Am 23.02.2016 um 22:22 schrieb Chris Reudenbach:

Augustin

just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are 
getting the information about the corner coordinates the subdatasets 
and so on. Together with Dominiks suggestion you can do something like 
this:


library(gdalUtils)
library(raster)
Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES")
wrffake<- "+proj=longlat +ellps=sphere +a=637 +b=637 +units=m 
+no_defs"

x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
  of="GTiff",
  ot="Byte",
  output_Raster=TRUE,
  verbose=TRUE)
wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619)
extent(x) <- extent(wrfExt)
projection(x) <- wrffake
plot(x)

Some remarks:
(1) I just took the first pair of  coordinates as derived from 
gdalinfo("geo_em.d01.nc")
you will find 4 different coordinate pairs (i did not proof which one 
is right


The data is staggered (as outlined by Dominik) So some of the corner 
coordinates belongs to the staggered data and the others coordinates  
to the unstaggered ones. You will find them marked


If you have installed the netcdf libs you easily can use ncview 
geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more 
information of the header.


Hope this helps

cheers Chris



Am 23.02.2016 um 21:11 schrieb Agus Camacho:

Thanks for that Dominik,

Giving that projection to either the locations, the raster layer 
generated

from the .nc file, or both, still did not work. I keep having locations
that should be on land falling far on the sea. Might this be a problem
derived from using raster with a file whose original grid distances 
are not

constant?

Here is a link with the original file which has the original coordinate
data.

https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0


2016-02-23 12:07 GMT-07:00 Dominik Schneider 
<dominik.schnei...@colorado.edu

:
This looks like WRF <http://www.wrf-model.org/index.php> data. I just
dealt with this.
The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
+a=637 +b=637 +units=m

+proj=lcc which is usually what wrf is run with.
The tricky part is:
+lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
because every WRF run is different (the WRF Preprocessing System 
optimizes

the projection for the domain).
and then there is probably no shift so you need(?) +x_0=0 +y_0=0

This gives:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 
+ellps=sphere

+a=637 +b=637 +units=m +no_defs

But, wrf users like to give out lat and  long so you need to assign it:
+proj=longlat +ellps=sphere +a=637 +b=637 +units=m +no_defs

and then reproject th

Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

2016-02-23 Thread Chris Reudenbach 

Augustin

just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are 
getting the information about the corner coordinates the subdatasets and 
so on. Together with Dominiks suggestion you can do something like this:


library(gdalUtils)
library(raster)
Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES")
wrffake<- "+proj=longlat +ellps=sphere +a=637 +b=637 +units=m 
+no_defs"

x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
  of="GTiff",
  ot="Byte",
  output_Raster=TRUE,
  verbose=TRUE)
wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619)
extent(x) <- extent(wrfExt)
projection(x) <- wrffake
plot(x)

Some remarks:
(1) I just took the first pair of  coordinates as derived from 
gdalinfo("geo_em.d01.nc")
you will find 4 different coordinate pairs (i did not proof which one is 
right


The data is staggered (as outlined by Dominik) So some of the corner 
coordinates belongs to the staggered data and the others coordinates  to 
the unstaggered ones. You will find them marked


If you have installed the netcdf libs you easily can use ncview 
geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more 
information of the header.


Hope this helps

cheers Chris



Am 23.02.2016 um 21:11 schrieb Agus Camacho:

Thanks for that Dominik,

Giving that projection to either the locations, the raster layer generated
from the .nc file, or both, still did not work. I keep having locations
that should be on land falling far on the sea. Might this be a problem
derived from using raster with a file whose original grid distances are not
constant?

Here is a link with the original file which has the original coordinate
data.

https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0


2016-02-23 12:07 GMT-07:00 Dominik Schneider  data. I just
dealt with this.
The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
+a=637 +b=637 +units=m

+proj=lcc which is usually what wrf is run with.
The tricky part is:
+lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
because every WRF run is different (the WRF Preprocessing System optimizes
the projection for the domain).
and then there is probably no shift so you need(?) +x_0=0 +y_0=0

This gives:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
+a=637 +b=637 +units=m +no_defs

But, wrf users like to give out lat and  long so you need to assign it:
+proj=longlat +ellps=sphere +a=637 +b=637 +units=m +no_defs

and then reproject the lat/long to lcc coordinates using this string:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
+a=637 +b=637 +units=m +no_defs

One word of caution, make sure you received the correct coordinates. Some
variables are run cell center while some are run at cell edge. It looks
like from your .nc file it was made by your collaborator so I assume they
are right.

That said, another word of caution, I found that the XLAT and XLONG
variables from WRF output aren't very precise. There is a "geogrid" file
from the preprocessing system that has the domain corners, resolution, nrow
and ncol from which you can make a better grid using the native projection
system (in my case it was a 4km grid). I suggest you try to get those.

I hope this helps... I have to run but wanted to save people too much head
scratching. I can get you running with more help tonight if you need.
Dominik


On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho 
wrote:


Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
more giant crab this time.

The creator of the .nc file also looked at this webpage:
http://www.pkrc.net/wrf-lambert.html
It seemed like the right proj4 string might be this one:

+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
 +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs

However this projection also does not allow me to adequately plot the
locations on the raster.

Here is the .nc file. it contains several layers.


https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0



2016-02-23 2:25 GMT-07:00 Michael Sumner :


On Tue, 23 Feb 2016 at 20:09 Roger Bivand  wrote:


On Tue, 23 Feb 2016, Alex Mandel wrote:


I made an attempt at it too. Investigating the original data, I'm

not

sure that the projection information supplied is correct for the

data

linked. When I load up the data in a unprojected space, the

coordinates

don't look at all similar to any Lambert projected data I have, they
actually look like Lat/Lon in some unprojected coordinate system,
perhaps a different spheroid than expected.

Does anyone have a link to the original data? Is is possible that

this is

the General Oblique Transformation used by modellers - that is

something

that feels like longlat but is

Re: [R-sig-Geo] calculate "regional" slope

2016-02-20 Thread Chris Reudenbach 
Dominik,

probably I did not catch exactly what  are your issues.

If you want to perform a multiscale analysis of the digital relief you 
should not focus on a 4,8,16 or more neighbors analysis. This would be 
helpful for a least cost path analysis or a single scale but multi 
neighborhood analysis.

If you have problems with your OS and just for getting an idea  you may 
use landserf  the great java tool of Jo Wood http://www.landserf.org/

As far as I understand your question I think a straightforward way would 
be to resample your elevation data in several steps to a coarser 
resolution and derive the information (aspect, watershed) you need. If 
you really want to have the information "on the other side of a 
mountain" You have to define what other side means. Most commonly it 
means something as behind the highest ridge.

But If you want to perform typical hydrological or morphometric analysis 
of the relief I would strongly suggest  the "mature" GIS software 
packages as GRASS7  SAGA GIS 2.12+ or also PCraster. All of them have 
very fast and very  well validated and reliable tools. PCRaster is more 
optimized with respect to the dynamics, while SAGA GIS and GRASS are 
more sophisticated in order to the implemented algorithms. GRASS is most 
optimized with respect to time consuming calculations. So if you "just" 
want to derive the information of the watersheds on different scales you 
can easily use all of these software packages.

There are some more aspects you should mention. If you need to analyze 
let's say 8 tiles of 5 deg SRTM data R is not able to deal with it.  
Even with PCRaster and SAGA you will need a lot of patience.

Maybe a final thought: we are dealing in students courses since some 
years using R and the mentioned GIS packages on all major OS platforms. 
Instead of struggling all the time with not satisfying work arounds for 
OSX you may think about a second partitions with Linux (or even 
Windows).  In the end this will save a lot of time...

If you need any support with one the mentioned topcis please let me know.

cheers Chris R


Am 20.02.2016 um 08:49 schrieb chris english:
>
> Dominik,
> There's mention of grass not working on osx 10.11 and a work around
> grassmac.wikidot.com 
>
> I was generally using PCRaster when I was doing watershed/flood 
> re-creations to handle the watershed delineations and the like and 
> don't have an opinion on Grass.
>
> Yeah, that transect was maybe a silly idea but it generalizes surface 
> to line, sort of like Mauana volcano is a cake you've cut in half. And 
> the line is outline.
>
> Chris
>
> On Feb 20, 2016 2:09 AM, "Dominik Schneider" 
>  > wrote:
>
> Chris R,
> I'd be happy to use grass7 but I can't get it to run on my mac.
> osx 10.11.3. I have a working qgis from kyngchaos, installed
> pandoc and cairo on top of that, and disabled sip but the
> grass7.app just hangs when I try to open it.
> Chris E,
> I will try to clarify. Rather than thinking about "side of a
> mountain", think about "side of a mountain **range**".  The point
> of calculating a regional slope is that if I am on the west side
> of a continental divide but on the east side of mountain, the
> function still returns a value indicating west-facing.  so maybe
> it's easier to think in meters.   lets assume my DEM is a 500m
> grid. the slope calculation would give a local value at 500m. This
> local slope might be east facing, but I am interested in the
> overall slope across 10km to indicate which way e.g. the watershed
> is draining.
> What do you mean with "transect along the z. I.e. roll your dem on
> it's side."?
>
>
> On Fri, Feb 19, 2016 at 2:50 PM, chris english
>  > wrote:
>
> Dominik,
> Sorry, I'm still trying to understand your 0.05 to 1.5 degree
> part of your problem.
>
> Otherwise, I think you are limited to 8 neighbors as this
> reflects the documentation as I read it.
>
> Even Roualt perhaps would be up in arms; but, there's nothing
> saying you can't do a 16 vs 8 neighbor. You'd have to examine
> the impacts thereafter, but basically you'd be amending some
> gdal (probably a line or two of code) for your purposes.
>
> There are a bunch of things to consider, theoretical and
> practical; but, why 16 better than 8. And more importantly as
> you relax, as a matter of rings (in this case), would your
> analytical result be better? Or potentially have any remainder
> meaning at all, I.e. I don't know my neighbor's neighbor's
> neighbor (does that get us out to 16?).
> And so generalizing beyond some given point might yield not
> much on a knn influence/likeness basis.
>
>

Re: [R-sig-Geo] calculate "regional" slope

2016-02-19 Thread Chris Reudenbach 
Dominik,

If you want to deal with bigger data sets in a fast and more flexible way I 
strongly suggest to use GRASS7. Rgrass7 provides a very good "r-ish" wrapper to 
use it from R. If You are interested i 'll post a typical setup and example for 
your question. Please note you need obligatory a GRASS 7 installation on your 
system.

Cheers Chris 


Am 19. Februar 2016 16:55:32 MEZ, schrieb Dominik Schneider 
:
>Thanks for the suggestion Chris.  I'm familiar with gdaldem, which
>raster::terrain is based on, to compute slope from a dem. I now realize
>that my example isnt a good one because neighbors=8 would achieve what
>I
>described. However I actually want some flexibility such that I can
>specifiy neighbors=16 so that it uses the next "ring" of cells.
>
>I played around with focal() with weight argument =
>matrix(rep(c(1,0,0,0,1),5),byrow=T) but couldn't figure out how to
>solve
>for a directional slope.
>
>On Fri, Feb 19, 2016 at 4:09 AM, chris english <
>englishchristoph...@gmail.com> wrote:
>
>> Dominik,
>>
>> r <- raster(nrows=22, ncols=20, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
>>  vals <- sample.int(1e3,440)
>> r[ ] <- vals
>> #raster::terrain
>> terr_r <- terrain(r, opt='slope', unit='degrees', neighbors=8)
>> Ah, but it appears you want up sampling to 1.5 degrees rather than
>0.5 deg.
>> so maybe spatial.tools::projectRaster_rigorous then raster:terrain.
>>
>> I'm inclined to end that last so maybe with a question mark. Sorry
>for an
>> essentially inconclusive response but I was happy to find terrain in
>any
>> case.
>> Chris
>>
>> On Fri, Feb 19, 2016 at 2:59 AM, Dominik Schneider <
>> dominik.schnei...@colorado.edu> wrote:
>>
>>> I need to calculate slope at different scales. In the case below, r
>is a
>>> 0.5deg resolution raster and I want the slope for 1.5 deg centered
>on each
>>> of those 0.5 deg pixels. I'm trying to estimate which side of
>mountain
>>> range each pixel is on. So the resulting raster would have the same
>number
>>> of pixels as r. The edges can be NA.
>>> any suggestions would be appreciated. Thanks
>>>
>>>
>>> r <- raster(nrows=22, ncols=20, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
>>> setValues(r,rnorm(440))
>>>
>>> [[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

-- 
Diese Nachricht wurde von meinem Android-Mobiltelefon mit K-9 Mail gesendet.
[[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] different projection transformation R and gdal commandline

2016-02-15 Thread Chris Reudenbach 

Hi Dominik,

If you use the gdalUtils package  there is no significant difference in 
the results using CLI or  R:


library(gdalUtils)
gdaltransform(s_srs="+proj=longlat +datum=WGS84", t_srs="+proj=lcc 
+lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0 +y_0=0  
+ellps=sphere  +a=637 +b=637 +units=m 
+no_defs",coords=matrix(c(-112.25,33.0), ncol = 2))

 [,1]  [,2] [,3]
[1,] -1306676 -629522.50


If I use  spTransform I can reproduce your results:


loc <- c(-112.25,33.0)
loc <- data.frame(matrix(unlist(loc), nrow=1, 
byrow=T),stringsAsFactors=FALSE)

colnames(loc)<-c("lon","lat")
coordinates(loc)<- ~lon+lat
proj4string(loc)<- CRS("+proj=longlat +datum=WGS84 +no_defs")
loc<-spTransform(loc,CRS('+proj=lcc +lat_1=28 +lat_2=50 
+lat_0=39.70001220694445 +lon_0=-98 +x_0=0 +y_0=0 +ellps=sphere 
+a=637 +b=637 +units=m +no_defs'))

loc@coords
   lon lat
1 -1306471 -629424

I suggest to focus on the sptransform() function

cheers Chris

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.3 LTS

locale:
 [1] LC_CTYPE=de_DE.UTF-8   LC_NUMERIC=C LC_TIME=de_DE.UTF-8
LC_COLLATE=de_DE.UTF-8
 [5] LC_MONETARY=de_DE.UTF-8LC_MESSAGES=de_DE.UTF-8 
LC_PAPER=de_DE.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C 
LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C


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

other attached packages:
[1] gdalUtils_2.0.1.7 raster_2.5-2  sp_1.2-2

loaded via a namespace (and not attached):
 [1] rgdal_1.1-3   tools_3.2.3   Rcpp_0.12.3 R.methodsS3_1.7.0 
codetools_0.2-14  grid_3.2.3
 [7] iterators_1.0.8   foreach_1.4.3 R.utils_2.2.0 
R.oo_1.19.0   lattice_0.20-33





Am 15.02.2016 um 21:37 schrieb Dominik Schneider:

Hi,
I'm struggling to use a custom projection. I am seeing differences with
someone using python proj4 bindings and when I compared my R results with
my commandline results I got even more confused. the coordinate
transformation is different for the two different methods.

could someone explain to me which one is wrong and why?
Thanks

R:
e=c(-112.25,-104.125,33,43.75)
box=as(extent(e),'SpatialPolygons')
proj4string(box)='+proj=longlat +datum=WGS84'
pstring='+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98
+x_0=0 +y_0=0 +ellps=sphere +a=637 +b=637 +units=m +no_defs'
xmin(spTransform(box,CRS(pstring)))
## [1] -1306471
ymin(spTransform(box,CRS(pstring)))
## [1] -713442.3

commandline:
iMac:~ $ gdaltransform -s_srs "+proj=longlat +datum=WGS84" -t_srs
"+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0
+y_0
=0  +ellps=sphere  +a=637 +b=637 +units=m +no_defs"
-112.25 33.
-1306675.75472246 -629522.472322824 0

[[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] LiDAR LAS file import into R via RSAGA

2016-02-12 Thread Chris Reudenbach 

Hi Clay,

RSAGA provides a comfortable platform independent system wrapper for 
SAGA but as stated  no two-way bindings.


If you want to make use of SAGA algorithms () it is IMHO the most 
straightforward way to import and export the raster files using the e.g. 
gdalUtils. If you use raster() directly please note that it may produce 
problems writing SAGA files directly. To avoid format complications  I 
use  most of the time GTiff. You can easily integrate the power of 
rLidar in this workaround concept.

A typical work flow may look like:


 # gdalwarp is used to convert the data format to SAGA
 # you may assign the projection information to the data (like example 
or change it with t_srs=)
  gdalwarp(file.in, file.out, overwrite=TRUE, 
s_srs=paste0('EPSG:',epsg.code), of='SAGA')


 # do something in SAGA using the files on the disk storage as "interface"
 # note SAGA works much more stable under differnet platforms when you 
provide the env settings


  # e.g. calculate wood's terrain indices   wood= 
1=planar,2=pit,3=channel,4=pass,5=ridge,6=peak

rsaga.geoprocessor('ta_morphometry',"Morphometric Features",env=myenv,
   list(DEM='mp_dem.sgrd',
FEATURES='mp_wood.sgrd',
SLOPE='mp_slope.sgrd',
LONGC='mp_longcurv.sgrd',
CROSC='mp_crosscurv.sgrd',
MINIC='mp_mincurv.sgrd',
MAXIC='mp_maxcurv.sgrd',
TOL_SLOPE=14.00
TOL_CURVE=0.1
EXPONENT=0.00
ZSCALE=1.00))

# read whatever you wand directly using raster()
peak.area<-raster('mp_wood.sdat')

For using SAGA and/or GRASS in a "mixed" environment with R you may find 
some help at our course pages:

 
http://moc.environmentalinformatics-marburg.de/doku.php?id=courses:msc:advanced-gis:lecture-notes:ag-ln-06

cheers Chris


Am 12.02.2016 um 07:51 schrieb Clay S:

Thanks Michael Sumner

rLidar works well.  But how do I pass the lidar data to another SAGA module
within R without saving the output of each step as a file in a whole series
of SAGA module calls from R?

I would like something like the dplyr %>% pipe command.

Thanks
Clay

On Fri, Feb 12, 2016 at 7:44 PM, Michael Sumner  wrote:


Try rLidar::readLAS as a straightforward alternative, it's not bullet
proof but works fine for a lot of examples. Happy to help extend it if your
format is not supported.

Cheers, Mike

On Fri, 12 Feb 2016 7:41 am Clay S  wrote:


Hi

Is it possible to import a Lidar LAS file into R with the
rsaga.geoprocessor function of the RSAGA package.  The following code runs
without an error, but all I get in R is the SAGA console output.

The R code below does save a SAGA point cloud file called "test" in my
working directory.  If I delete the `POINTS = "test"` parameter, then
nothing happens as expected.  However, I want to work with the data in R,
and pass the point cloud to other SAGA modules. This can be done in SAGA
without having to save the point cloud.  Can the same be done from R?

How can I get the point cloud imported into the R environment?  Do I need
to coerce it into a matrix, data frame, or Spatial Points object?


R code:

###

library(RSAGA)


# Set RSAGA environment.

env <- rsaga.env()


# SET WORKING DIRECTORY FIRST!!

dir <- getwd()


# Get module parameters

lasfilelist <- list.files(dir, pattern =".las$", full.names=TRUE,
recursive=FALSE)


# send command to RSAGA

test <- rsaga.geoprocessor(lib = "io_shapes_las", module = "Import LAS
Files",
param = list(FILES = lasfilelist[1], POINTS = "test",
 T=TRUE, i=TRUE, a=FALSE, r=TRUE, c=FALSE,
u=FALSE,
 n=TRUE, R=FALSE, G=FALSE, B=FALSE,
e=FALSE,
d=FALSE,
 p=FALSE, C=FALSE, VALID=FALSE,
RGB_RANGE=FALSE),
env = env)

###


Here is the console output (the test variable stores the same thing):

###

library path: C:\OSGEO4~1\apps\saga\modules\io_shapes_las.dll
library name: Import/Export - LAS
tool name   : Import LAS Files
author  : O. Conrad, V. Wichmann (c) 2009

Parameters

Input Files: "C:/Users/whcms0/'.CPT/Tarewa/Lidar/AX300711.las"
Point Clouds: No objects
gps-time: yes
intensity: yes
scan angle: no
number of the return: yes
classification: no
user data: no
number of returns of given pulse: yes
red channel color: no
green channel color: no
blue channel color: no
edge of flight line flag: no
direction of scan flag: no
point source ID: no
rgb color: no
Check Point Validity: no
R,G,B value range: 16 bit

Parsing AX300711.las ...
Save point cloud: test...




Here is my env:




env

$workspace

[1] "."


$cmd

[1] "saga_cmd.exe"


$path

[1]

Re: [R-sig-Geo] raster pkg: overwriting the dimensions detected by netcdf4 (PERSIANN-CDR)

2016-02-03 Thread Chris Reudenbach 
If you look at the file with ncview and analyse the header with ncdump - 
h you will find that en extent is missing.

The standard tool ncview shows also a "flipped version"

Try to use gdal more directly like this snippet

library(raster)
library(gdalUtils)
library(mapview)

# necessary because of nc file see configuration options 
http://www.gdal.org/frmt_netcdf.html

Sys.setenv(GDAL_NETCDF_BOTTOMUP="NO")
extentNC <- extent(0, 360,-60, 60)
proj4str <- '+proj=longlat +datum=WGS84 +no_defs'
ncfilename <- "PERSIANN-CDR_v01r01_20090101_c20140523.nc"
dataset<- "precipitation"
band <-1

x<-gdal_translate(paste0("NETCDF:",ncfilename,":",dataset), 
paste0(ncfilename,".tif"),

   b=band,
   of="GTiff",
   output_Raster=TRUE,
   verbose=TRUE,
   overwrite=TRUE)
# due to GDAL_NETCDF_BOTTOMUP="NO" we have to flip the raster
x <- flip(x, direction='y')
extent(x) <- extentNC
projection(x) <- proj4str

mapview(x)
plot(x)


cheers Chris


Am 03.02.2016 um 16:20 schrieb MAURICIO ZAMBRANO BIGIARINI:

On 3 February 2016 at 11:34, Michael Sumner  wrote:


On Thu, 4 Feb 2016 at 01:26 MAURICIO ZAMBRANO BIGIARINI
 wrote:

Dear spatial community,

While reading some netCDFfiles with the raster package , I got a
"rotated" file, where the columns where read as rows and vice versa:

- START 
x <- raster("PERSIANN-CDR_v01r01_20090101_c20140523.nc")
Loading required namespace: ncdf4
plot(x)
x
class   : RasterLayer
dimensions  : 1440, 480, 691200  (nrow, ncol, ncell)
resolution  : 0.25, 0.25  (x, y)
extent  : -60, 60, 0, 360  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
names   : NOAA.Climate.Data.Record.of.PERSIANN.CDR.daily.precipitation
z-value : 2009-01-01
zvar: precipitation
- END 


Transpose might be sufficient?

tx <- t(x)

plot(tx)
tx


Thanks Michale,

In fact, that is what I'm doing for reading the files, and the spatial
extent is correct and the precipitation values also seems to be right.

However, my question was more "conceptual", in the sense that the
provider says that the original file is NOT rotated, while the only
way  I have to get the correct files is rotating the original file
with the 't' command you mentioned.

In addition, while reading the file with the latest version of GRASS
GIS (7.0.3), and I got the following error message:

"
(Tue Feb  2 19:41:29 2016)
r.in.gdal input=/home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
output=PERSIANN_CDR_v01r01_20090101_c20140523 -e
Warning 1: dimension #2 (lat) is not a Longitude/X
dimension.
Warning 1: dimension #1 (lon) is not a Latitude/Y dimension.
ERROR: Input raster map is flipped or rotated - cannot import. You may
use 'gdalwarp' to transform the map to North-up.
"


So, I'm quite sure that the original file IS rotated, but I didn't
have more technical arguments to give to the data provider to
demonstrate that the file is rotated, because they insist that the
rotation is because R is not reading the file in the proper way

mzb

=
"In the end, it's not the years in your life that count.
  It's the life in your years". (Abraham Lincoln)
=
Linux user #454569 -- Linux Mint user


60W-60E and 0-360N, while the correct extent should be: 60S, 60N, 0E,
360E.

The file can be downloaded from:

ftp://data.ncdc.noaa.gov/cdr/persiann/files/2009/PERSIANN-CDR_v01r01_20090101_c20140523.nc


After contacting the providers of the file, they mentioned that when
reading the file I have to provide the dimension information, i.e.,
480 rows and 1440 columns, while the raster package detects the other
way around.

I would highly appreciate if you could tell me:

1) is the previous situation a bug of my netCDF driver (netCDF 4.1.3,
GDAL 1.11.2, released 2015/02/10) or a problem of the dimensions
actually recorded in the netCDF file ?

2) is there any way to  overwrite the dimensions detected by netcdf4
when reading the file with the 'raster' command ?


My sessionInfo():
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.2 LTS

locale:
  [1] LC_CTYPE=en_GB.UTF-8   LC_NUMERIC=C
  [3] LC_TIME=en_GB.UTF-8LC_COLLATE=en_GB.UTF-8
  [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_GB.UTF-8
  [7] LC_PAPER=en_US.UTF-8   LC_NAME=C
  [9] LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] raster_2.5-2 sp_1.2-1

loaded via a namespace (and not attached):
[1] tools_3.2.3 Rcpp_0.12.3 ncdf4_1.15  grid_3.2.3
[5] lattice_0.20-33


Thanks in advance and kind regards

Mauricio Zambrano-Bigiarini, PhD