Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 2 hours 32 min ago

ApproxNA producing different and incorrect results in stacks of large rasters

Thu, 03/02/2017 - 07:52
Hello, I’m trying to confirm if there is a bug or at least unexplained
change in behavior in the function approxNA in the most recent version of
the R Raster package (Raster 2.5-8).



I have a small reproducible example below. The original symptoms (too large
to share) were that most the output of the approxNA command appeared to be
the same as the input – it had lots and lots of NA values that had not been
interpolated… but when I switched back to an older version of the package
(2.3-12) and an older version of R (3.3.1) it returns the correct
interpolated results.  I have tried it with raster datatype set to INT2S
and with no particular raster datatype set.



I have checked the documentation in case there is some reason that I should
expect the results to differ but there are no changes between the old
package version that worked correctly and the new package version.



In my reproducible example there are fewer errors but they are quite
obvious to see just from the plotting view (included in the example). In
these plots you can see that the northern parts of several layers have NA
values remaining that should have been interpolated.



Official system bug report details are shown at the bottom of this email;
I’m running 64-bit Windows 7 and 64-bit R 3.3.2 with raster 2.5-8 with an
R-Studio interface, but I have tested also on R-console with the same
results and on a different computer that has the same configuration, I have
also tested in 32-bit R with the same results.



A trigger seems to be the presence of NA values in the first or last layer,
but the errors occur in cells that do not have NA values in the first or
last layer: for example the top left cells are only missing in layer 2 in
the original stack, so it should be a very simple interpolation for those
cells from layer 1 to layer 3.



I did follow the R instructions for submitting a bug report a couple of
weeks ago but I’m not sure if that contact info is correct, or if it might
have gone into a spam folder.



Best regards, and thanks for any suggestions



Jodi



##### start  Reproducible example



library(raster)

# make a large raster, it must be quite big to have a problem

r <- raster(ncols=8142, nrows=3285)

# make six rasters to stack, use different ranges to better see the
interpolation results

r1 <- setValues(r, round(10 * runif(ncell(r)),0))

r2 <- setValues(r, NA)

r3 <- setValues(r, round(100 * runif(ncell(r)),0))

r4 <- setValues(r, round(1000 * runif(ncell(r)),0))

r5 <- setValues(r, round(50 * runif(ncell(r)),0))

r6 <- setValues(r, round(100 * runif(ncell(r)),0))



# insert some NA values

r1[100:600,] <- NA

r3[,1500:1950] <- NA

r5[,400:800] <- NA

r6[2500:2900,] <- NA



# insert some constant values to make it easier to see if interpolation is
working

r4[300:500,] <- 750

r6[300:400,] <- 20

r6[400:500,] <- 100

# make a stack

s <- stack(r1,r2,r3,r4,r5,r6)



# see what the pre-interpolation stack looks like

plot(s)



# interpolate, this takes a while

x2 <- approxNA(s, method = "linear", rule=2)



# see what the post interpolation looks like, there should be filled in
values for all parts of all layers

# but instead it’s retaining many of the NA values

plot(x2)



## End reproducible example







--please do not edit the information below--

Package: raster
 Version: 2.5-8
 Maintainer: Robert J. Hijmans <[hidden email]>
 Built: R 3.3.1; x86_64-w64-mingw32; 2016-07-19 01:23:33 UTC; windows

R Version:
 platform = x86_64-w64-mingw32
 arch = x86_64
 os = mingw32
 system = x86_64, mingw32
 status =
 major = 3
 minor = 3.2
 year = 2016
 month = 10
 day = 31
 svn rev = 71607
 language = R
 version.string = R version 3.3.2 (2016-10-31)
 nickname = Sincere Pumpkin Patch

Windows 7 x64 (build 7601) Service Pack 1

Locale:
 LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

Search Path:
 .GlobalEnv, package:raster, package:sp, package:stats,
 package:graphics, package:grDevices, package:utils, package:datasets,
 package:methods, Autoloads, package:base







[hidden email]
Landscape and Quantitative Ecologist
Inventory and Monitoring Program
Southern Colorado Plateau Network
(928)-523-6142


National Park Service SCPN
Northern Arizona University
525 S Beaver, Bldg 20, Rm 131
PO Box 5663
Flagstaff, AZ, 86011-5663

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [FORGED] spsample: defining the size of the cell in an hexagonal grid

Thu, 03/02/2017 - 07:01
Thank you very much Rolf.

Manuel

2017-03-01 2:03 GMT-06:00 Rolf Turner <[hidden email]>:

> On 01/03/17 17:33, Manuel Spínola wrote:
>
>> Dear list members,
>>
>> How can I define the size in hectares of the cell when doingan hexagonal
>> grid with the function spsample.  What it means "cellsize"?  I want to
>> create cell sizes according to some specific number of hectares.  How can
>> I
>> do it?
>>
>>
>> data(meuse.grid)
>> gridded(meuse.grid) = ~x+y
>> plot(meuse.grid)
>> HexPts <-spsample(meuse.grid, type="hexagonal", cellsize=1000)
>> HexPols <- HexPoints2SpatialPolygons(HexPts)
>> plot(HexPols[meuse.grid,], add=TRUE)
>>
>
>
> I have no knowledge of spsample() nor of its argument cellsize (doesn't
> the help for spsample() tell you what this argument means?) but it might be
> useful for you to realise (if you are not already aware of this) that the
> area of a regular hexagon is
>
>     A = 3*sqrt(3)*s^2/2
>
> where s is the length of side of the hexagon.  So if you want hexagons of
> area 42 hectares you would take their sides to be of length 402.0673 metres.
>
> cheers,
>
> Rolf Turner
>
> --
> Technical Editor ANZJS
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>


--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email] <[hidden email]>
[hidden email]
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Suggestions for irregular grid netCDF data?

Wed, 03/01/2017 - 22:43
Hi there,

my general advice is

* use raster to read the data layers (including the coordinates, which in
this context are just data)
* treat each layer "index space", i.e. setExtent(r, extent(0, nrow(x), 0,
ncol(x)) (or possibly transpose for some NetCDF files)
* use raster's cell-index logic to do extractions and aggregations in that
index space, it's neat because cropping inherits the parent georeference so
it's "global"
* use resampling to put real-world data objects into that index space for
extractions (you have to worry about the resolution and the boundary to do
that)

I have a few R packages that do this for ROMS and CMIP5 that can be
generally applied, but there not in great shape atm.

This rough document shows some of what you can do, it was written for a
particular purpose I had exploring leaflet which is on the backburner for
now, so it doesn't really go anywhere.

http://rpubs.com/cyclemumner/roms0

I'm happy to help you get what you need though, we deal with scads of this
kind of data and it doesn't fit any good mould except the "code it
yourself" club. It's not a huge amount of code but the "many coordinate
systems" can be bewildering and it's really specifically bound to the
intricacies of raster and sp and rgdal and ncdf4.

I'll try to look at your examples sometime soon.

Cheers, Mike.


On Thu, 2 Mar 2017 at 14:53 Michael Treglia <[hidden email]> wrote:

> ​Hi all,
>
> I've been trying to work with a somewhat complex netCDF dataset, and wanted
> to see if anybody has suggestions.
>
> The data are for marine variables in the form of a remote netCDF served
> through OPENDAP. One of the complexities is that it comes in an irregular
> grid, as fine as 25m and as coarse as 400m (coarser as you move further
> from the coast) - this variable spacing makes it tough to work with in the
> raster package, which seems like it would be the easiet.  Is there a
> straightforward, relatively simple way to work with a dataset like this,
> getting entire layers at a time? The data have a slew of different
> variables, some of which also have time steps and different depth horizons.
>
> So far, I’ve been able to connect to the data using the ncdf4 package in R,
> and load simple, 3 dimensional data (e.g., depth). I can also force layers
> to load using the raster package, adding the argument
> ‘stopIfNotEqualSpaced=FALSE', though that obviously has issues with the
> irregular grid (though if it worked, would make it easier to pull
> information by date and depth horizons). I'm having a tough time getting
> data via the ncdf4 package for specific dates/depth horizons though. [note:
> I end up having to use R in a linux environment for this, as I understand
> the windows implementation for is more limited, and throws an ‘Unknown file
> format’ error.]
>
> The data are here:
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
> ,
> and here’s the web-interface for the data:
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc.html
> .
> I've included some sample R code below - and am totally open to alternative
> ways that might be better for working with these data [still new-ish to
> netCDFs, and have mainly dealt with much simpler, soil or climate
> datasets].
>
> Thanks!
> Mike
>
>
> Sample R Code:
> #Load and work w/ data in ncdf4
> library(ncdf4)
>
> nyhops.nc <- nc_open("
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
> ")
>
> #Get Point Data with Depth Measurements
> depth.nc <- ncvar_get(nyhops.nc, 'depth')
> lat.nc <- ncvar_get(nyhops.nc, 'lat')
> long.nc <- ncvar_get(nyhops.nc, 'lon')
>
> depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
> depth=as.vector(depth.nc)) #Can then be exported, or converted to spatial
> object
>
> #Trying to get the entirety of a single time point and depth horizon, but
> not quite getting the start/count args - this gets a single data point it
> seems
> temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(10,10,1,1),
> count=c(1,1,1,1))
>
>
>
> #Load layers w/ raster package.
> library(raster)
> nyhops.url <- "
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
> "
> depth <- raster(nyhops.url, varname="depth", lvar=1, level=1,
> stopIfNotEqualSpaced=FALSE) #Forcing it to read layer despite unequal
> spacing of coords
> plot(depth) #Plot looks distorted, as expected given irregular grid
> spacing. [I realize the data aren't projected at this point]
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Suggestions for irregular grid netCDF data?

Wed, 03/01/2017 - 21:53
​Hi all,

I've been trying to work with a somewhat complex netCDF dataset, and wanted
to see if anybody has suggestions.

The data are for marine variables in the form of a remote netCDF served
through OPENDAP. One of the complexities is that it comes in an irregular
grid, as fine as 25m and as coarse as 400m (coarser as you move further
from the coast) - this variable spacing makes it tough to work with in the
raster package, which seems like it would be the easiet.  Is there a
straightforward, relatively simple way to work with a dataset like this,
getting entire layers at a time? The data have a slew of different
variables, some of which also have time steps and different depth horizons.

So far, I’ve been able to connect to the data using the ncdf4 package in R,
and load simple, 3 dimensional data (e.g., depth). I can also force layers
to load using the raster package, adding the argument
‘stopIfNotEqualSpaced=FALSE', though that obviously has issues with the
irregular grid (though if it worked, would make it easier to pull
information by date and depth horizons). I'm having a tough time getting
data via the ncdf4 package for specific dates/depth horizons though. [note:
I end up having to use R in a linux environment for this, as I understand
the windows implementation for is more limited, and throws an ‘Unknown file
format’ error.]

The data are here:
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc,
and here’s the web-interface for the data:
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc.html.
I've included some sample R code below - and am totally open to alternative
ways that might be better for working with these data [still new-ish to
netCDFs, and have mainly dealt with much simpler, soil or climate datasets].

Thanks!
Mike


Sample R Code:
#Load and work w/ data in ncdf4
library(ncdf4)

nyhops.nc <- nc_open("
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
")

#Get Point Data with Depth Measurements
depth.nc <- ncvar_get(nyhops.nc, 'depth')
lat.nc <- ncvar_get(nyhops.nc, 'lat')
long.nc <- ncvar_get(nyhops.nc, 'lon')

depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
depth=as.vector(depth.nc)) #Can then be exported, or converted to spatial
object

#Trying to get the entirety of a single time point and depth horizon, but
not quite getting the start/count args - this gets a single data point it
seems
temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(10,10,1,1),
count=c(1,1,1,1))



#Load layers w/ raster package.
library(raster)
nyhops.url <- "
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
"
depth <- raster(nyhops.url, varname="depth", lvar=1, level=1,
stopIfNotEqualSpaced=FALSE) #Forcing it to read layer despite unequal
spacing of coords
plot(depth) #Plot looks distorted, as expected given irregular grid
spacing. [I realize the data aren't projected at this point]

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

SpatialDataFrame List Union

Wed, 03/01/2017 - 12:05
Dear all,

I'm working on a code and I don't understand something probably related to
the writting rules of R

My code :

library(rgdal)

shp_paths = list.files(path = "X:/IGN/BDTOPO_V2.1", pattern
="TRONCON_COURS_EAU.SHP",
                       full.names = T,recursive = T,include.dirs = F)

shp_objects <- lapply(shp_paths, function(x) {readOGR(dsn=x,
                       layer=ogrListLayers(x))}) ### Liste de
SpatialLineDataFrames


"shp_objects" is a large list of 16 elements (192 Mb)

What I want to do is an union of these 16 elements.

Does anybody has a pist ?

Thanks !

Cheers

Tristan
Hydromorphologist Engineer
Seine Normandie Water Agency
Rouen, France
--
Tristan Bourgeois

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

EGU: R’s deliberate role in Earth sciences

Wed, 03/01/2017 - 10:59
This year's EGU meeting has a session with 16 presentation on "R’s
deliberate role in Earth sciences"; for the program see:

http://meetingorganizer.copernicus.org/EGU2017/pico/24971
--
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
signature.asc (484 bytes) Download Attachment

Re: split SLDF with attributes

Wed, 03/01/2017 - 08:15
*Hi Melanie,*

*Sorry for the delayed response. Your loop made exactly what I want it.
Thanks Mel, it's really helpful.*
*Now, I'm trying to manipulated the function last.merge to join my segments
by length. My aim is dividing the transects by 10km (n segments). However,
If the last segment is shorter than 5km, I have two options: *


*option 1)W*e can distribute the remainder among the n segments
(segments>10km).
raw-lines: 10+10+2
output-1:   11+11


*option 2)*  We can add 1 segment(n+1), and divided by segments smaller
than 10km.
raw-lines: 10+10+7
process: 27/3
output-2:     9+9+9

*I think my main aim is still far away. Then I decided to start with
something simple, the function **SegmentSpatialLines split the segments *
*effectively* *in the 10km and merge the last segment with the penultimate.
**I thought it would be easy define a cut-off to decide if the shorter
segments can be merge, but I was wrong.*

*1) First,* the transects were divided by length(10km) into segments.
segments:
10+10+6+10+3

*2) Second*, the length of transects is not a multiple of given
length(10km), if last segment is shorter we can merge it with penultimate
segment, if we wish(merge.last=T).
10+16+13

*3) Third*, I would like merge the sorter last segment with the penultimate
segment, only if the length is <5km.
segments:
10+10+null+13

*#I try to include a cut-off with these modifications:*
#*A)*      if ((merge.last && length(segments) <5000))#5km cut-off
##But the output of A was the same as if ((merge.last && length(segments)
>1)) and merge.last=F
#
I thought may be the length(segments) here is only to give you the nrow and
not the length of the segments. And I decided to use other ways to
calculate the length. The funcion gLenght of the rgeos package, it's not
possible because the segments in the function SegmentSpatialLines are in a
list. I also projected the segments on the loop before the last.merge, but
without success.  Then I try an other option with the total_lenght.
#*B)*
##change segmentspatilines function
SegmentSpatialLines <- function(sl, length = 0, n.parts = 0, merge.last =
F) {
 stopifnot((length > 0 || n.parts > 0))
 id <- 0
 newlines <- list()
 sl <- as(sl, "SpatialLines")
 for (lines in sl@lines) {
   for (line in lines@Lines) {
     crds <- line@coords
     # create segments
     segments <- CreateSegments(coords = crds, length, n.parts)
     # calculate total length line
       total_length <- 0
     for (i in 1:(nrow(crds) - 1)) {
       d <- sqrt((crds[i, 1] - crds[i + 1, 1])^2 + (crds[i, 2] - crds[i +
1, 2])^2)
       total_length <- total_length + d
     }

     if ((merge.last && length(segments) >1) &
(total_length(segments)<5000))

       {
       # in case there is only one segment, merging would result into error
       segments <- MergeLast(segments)
     }
     # transform segments to lineslist for SpatialLines object
     for (segment in segments) {
       newlines <- c(newlines, Lines(list(Line(unlist(segment))), ID =
as.character(id)))
       id <- id + 1
     }
   }
 }
 return(SpatialLines(newlines))
}




##But the output of B was the same as if ((merge.last && length(segments)
>1)) and merge.last=F

Thanks in advance,

Marta

>
>>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [FORGED] spsample: defining the size of the cell in an hexagonal grid

Wed, 03/01/2017 - 02:03
On 01/03/17 17:33, Manuel Spínola wrote:
> Dear list members,
>
> How can I define the size in hectares of the cell when doingan hexagonal
> grid with the function spsample.  What it means "cellsize"?  I want to
> create cell sizes according to some specific number of hectares.  How can I
> do it?
>
>
> data(meuse.grid)
> gridded(meuse.grid) = ~x+y
> plot(meuse.grid)
> HexPts <-spsample(meuse.grid, type="hexagonal", cellsize=1000)
> HexPols <- HexPoints2SpatialPolygons(HexPts)
> plot(HexPols[meuse.grid,], add=TRUE)

I have no knowledge of spsample() nor of its argument cellsize (doesn't
the help for spsample() tell you what this argument means?) but it might
be useful for you to realise (if you are not already aware of this) that
the area of a regular hexagon is

     A = 3*sqrt(3)*s^2/2

where s is the length of side of the hexagon.  So if you want hexagons
of area 42 hectares you would take their sides to be of length 402.0673
metres.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

spsample: defining the size of the cell in an hexagonal grid

Tue, 02/28/2017 - 22:33
Dear list members,

How can I define the size in hectares of the cell when doingan hexagonal
grid with the function spsample.  What it means "cellsize"?  I want to
create cell sizes according to some specific number of hectares.  How can I
do it?


data(meuse.grid)
gridded(meuse.grid) = ~x+y
plot(meuse.grid)
HexPts <-spsample(meuse.grid, type="hexagonal", cellsize=1000)
HexPols <- HexPoints2SpatialPolygons(HexPts)
plot(HexPols[meuse.grid,], add=TRUE)


--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email] <[hidden email]>
[hidden email]
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: rasterTopolygons: singleparts often needed

Sun, 02/26/2017 - 13:39
Yes, this is correct. I refer to integrating disaggregate as an option, as
I think that this is what the user needs in most cases.
Agus

On Sun, Feb 26, 2017 at 1:14 PM, John Baumgartner <[hidden email]> wrote:
> Not at my computer, but does sp::disaggregate(v) return your singlepart
> polygons?
>
> On Sun, 26 Feb 2017 at 22:40, Agustin Lobo <[hidden email]> wrote:
>
>> Could raster::rasterTopolygons() be modified in future versions
>>
>> to optionally output singlepart polygons? e.g, in the following
>>
>> r <-raster(nrow=16, ncol=16)
>>
>> d <- rep(rep(1:4,rep(4,4)),8)
>>
>> r[] <- c(d,rev(d))
>>
>> plot(r)
>>
>> v <- rasterToPolygons(r,dissolve=TRUE)
>>
>> plot(v[1,])
>>
>> dim(v@data)
>>
>>
>>
>> the new option
>>
>> v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)
>>
>> would result on
>>
>> plot(v[1,])
>>
>> displaying one single square and
>>
>> dim(v@data)
>>
>> would be 8 rows and 1 column
>>
>>
>>
>> (equivalent to gdal_polygonize)
>>
>>
>>
>> Agus
>>
>>
>>
>> _______________________________________________
>>
>> R-sig-Geo mailing list
>>
>> [hidden email]
>>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to draw a legend when plotting sf (simple feature) objects

Sun, 02/26/2017 - 12:17
Thank you very much Edzer,

it works for me:

> ggplot(nc) + geom_sf(aes(fill = SID79))
Warning message:
st_crs<- : replacing crs does not reproject data; use st_transform for that

2017-02-26 11:46 GMT-06:00 Edzer Pebesma <[hidden email]>:

>
>
> On 26/02/17 13:23, Manuel Spínola wrote:
> > Dear list members,
> >
> > How can I draw a legend when plotting an sf (simple feature) object.
> >
> > nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)
> >
> > plot(nc["SID79"])
> >
>
> You're in base plot, so at this moment you can add elements
> incrementally to the plot, see e.g. ?legend.
>
> I have no plans of making this easy or automatic, because I don't think
> that base plot is the right place for it (although packages raster,
> spatstat and fields try hard), but would happily consider contributions.
>
> geom_sf in ggplot2 (now in sf branch on github) will follow up
> sp::spplot for sf objects; see e.g.
>
> https://github.com/edzer/sfr/issues/88#issuecomment-276738460
>
> nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet =
> TRUE)
> ggplot(nc) + geom_sf(aes(fill = SID79))
>
> Trying this with the gpkg you read in above fails for me; it looks like
> ggplot now wrongly assumes that geometry columns are always called
> `geometry' (despite the docs), which is not the case.
> --
> Edzer Pebesma
> Institute for Geoinformatics  (ifgi),  University of Münster
> Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
> Journal of Statistical Software:   http://www.jstatsoft.org/
> Computers & Geosciences:   http://elsevier.com/locate/cageo/
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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
[hidden email] <[hidden email]>
[hidden email]
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: How to draw a legend when plotting sf (simple feature) objects

Sun, 02/26/2017 - 11:46


On 26/02/17 13:23, Manuel Spínola wrote:
> Dear list members,
>
> How can I draw a legend when plotting an sf (simple feature) object.
>
> nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)
>
> plot(nc["SID79"])
>

You're in base plot, so at this moment you can add elements
incrementally to the plot, see e.g. ?legend.

I have no plans of making this easy or automatic, because I don't think
that base plot is the right place for it (although packages raster,
spatstat and fields try hard), but would happily consider contributions.

geom_sf in ggplot2 (now in sf branch on github) will follow up
sp::spplot for sf objects; see e.g.

https://github.com/edzer/sfr/issues/88#issuecomment-276738460

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
ggplot(nc) + geom_sf(aes(fill = SID79))

Trying this with the gpkg you read in above fails for me; it looks like
ggplot now wrongly assumes that geometry columns are always called
`geometry' (despite the docs), which is not the case.
--
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
signature.asc (484 bytes) Download Attachment

How to draw a legend when plotting sf (simple feature) objects

Sun, 02/26/2017 - 06:23
Dear list members,

How can I draw a legend when plotting an sf (simple feature) object.

nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)

plot(nc["SID79"])



Best,

--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email] <[hidden email]>
[hidden email]
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: rasterTopolygons: singleparts often needed

Sun, 02/26/2017 - 06:14
Not at my computer, but does sp::disaggregate(v) return your singlepart
polygons?

On Sun, 26 Feb 2017 at 22:40, Agustin Lobo <[hidden email]> wrote:

> Could raster::rasterTopolygons() be modified in future versions
>
> to optionally output singlepart polygons? e.g, in the following
>
> r <-raster(nrow=16, ncol=16)
>
> d <- rep(rep(1:4,rep(4,4)),8)
>
> r[] <- c(d,rev(d))
>
> plot(r)
>
> v <- rasterToPolygons(r,dissolve=TRUE)
>
> plot(v[1,])
>
> dim(v@data)
>
>
>
> the new option
>
> v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)
>
> would result on
>
> plot(v[1,])
>
> displaying one single square and
>
> dim(v@data)
>
> would be 8 rows and 1 column
>
>
>
> (equivalent to gdal_polygonize)
>
>
>
> Agus
>
>
>
> _______________________________________________
>
> R-sig-Geo mailing list
>
> [hidden email]
>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Assign a unique ID number to each patch in a raster

Sun, 02/26/2017 - 06:09
Yes, my gist mentioned by Agustin should do the trick. See the example in
the comment at https://gist.github.com/johnbaums/6d155cfca28e02b05ad5.

On Sun, 26 Feb 2017 at 21:42, Agustin Lobo <[hidden email]> wrote:

> I think Nelly points out that package raster lacks a more general
>
> clumping tool (i.e., the equivalent to r.clump in grass) that I've
>
> often missed too.
>
> I'm not sure circumventing the way Ben suggest would work for large
>
> raster layers, for which package raster was designed.
>
> I have a note to myself regarding this
>
> https://gist.github.com/johnbaums/6d155cfca28e02b05ad5
>
> and
>
> www.openforis.org/tools/geospatial-toolkit.html
>
> but have not tried it.
>
>
>
> Agus
>
>
>
> On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper <[hidden email]> wrote:
>
> > Hi,
>
> >
>
> > I think you are asking for connected component labeling. raster::clump()
> expects your features to be like islands surrounded by NA or 0 background.
> There might be more Raster* centric ways to do it, but you could use plain
> old vanilla image processing using the very good imager package (
> https://cran.r-project.org/web/packages/imager/ <
> https://cran.r-project.org/web/packages/imager/> ).  The steps below show
> in color your actual values with the connected component IDs labeled on
> each cell.  I used the 8-connected labeling but 4-connected is also
> available.
>
> >
>
> > library(raster)
>
> > library(imager)
>
> >
>
> > set.seed(10)
>
> > nc <- 8
>
> > nr <- 8
>
> > MIN <- 1
>
> > MAX <- 4
>
> > LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")
>
> >
>
> > r <- raster(ncols=nc, nrows=nr)
>
> > r[] <- round(runif(ncell(r),MIN, MAX),digits=0)
>
> > cc <- r
>
> >
>
> > img <- imager::as.cimg(as.matrix(r))
>
> > labeled <- imager::label(img, high_connectivity = TRUE)
>
> > cc[] <- as.matrix(labeled)
>
> >
>
> > plot(r, col = LUT)
>
> > txt <- cc[]
>
> > xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))
>
> > text(xy[,1], xy[,2], txt )
>
> >
>
> >
>
> > Cheers,
>
> > Ben
>
> >> On Jan 26, 2017, at 4:48 PM, Nelly Reduan <[hidden email]> wrote:
>
> >>
>
> >> Hello,
>
> >> I would like to assign a unique ID number to each patch (a patch is
> composed of a set of adjacent cells) as shown in this figure:
>
> >>  <2328.png>
>
> >>
>
> >> I tested the clump function (package raster) by applying an eight
> adjacent cells rule but this function mixes all cell values into single
> patch.
>
> >>
>
> >> Here is a code example to create a raster. The raster contains values
> ranged from 1 to 9.
>
> >> r <- raster(ncols=12, nrows=12)
>
> >> r[] <- round(runif(ncell(r),1,9),digits=0)
>
> >> plot(r)
>
> >>
>
> >> Is there a way to assign a unique ID number to grouping cells that form
> a patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?
>
> >>
>
> >> Thanks a lot for your time
>
> >> Nell
>
> >>
>
> >>
>
> >> _______________________________________________
>
> >> R-sig-Geo mailing list
>
> >> [hidden email] <mailto:[hidden email]>
>
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>
> > Ben Tupper
>
> > Bigelow Laboratory for Ocean Sciences
>
> > 60 Bigelow Drive, P.O. Box 380
>
> > East Boothbay, Maine 04544
>
> > http://www.bigelow.org
>
> >
>
> >
>
> >
>
> >
>
> >         [[alternative HTML version deleted]]
>
> >
>
> > _______________________________________________
>
> > R-sig-Geo mailing list
>
> > [hidden email]
>
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> _______________________________________________
>
> R-sig-Geo mailing list
>
> [hidden email]
>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

rasterTopolygons: singleparts often needed

Sun, 02/26/2017 - 05:15
Could raster::rasterTopolygons() be modified in future versions
to optionally output singlepart polygons? e.g, in the following
r <-raster(nrow=16, ncol=16)
d <- rep(rep(1:4,rep(4,4)),8)
r[] <- c(d,rev(d))
plot(r)
v <- rasterToPolygons(r,dissolve=TRUE)
plot(v[1,])
dim(v@data)

the new option
v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)
would result on
plot(v[1,])
displaying one single square and
dim(v@data)
would be 8 rows and 1 column

(equivalent to gdal_polygonize)

Agus

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Assign a unique ID number to each patch in a raster

Sun, 02/26/2017 - 04:41
I think Nelly points out that package raster lacks a more general
clumping tool (i.e., the equivalent to r.clump in grass) that I've
often missed too.
I'm not sure circumventing the way Ben suggest would work for large
raster layers, for which package raster was designed.
I have a note to myself regarding this
https://gist.github.com/johnbaums/6d155cfca28e02b05ad5
and
www.openforis.org/tools/geospatial-toolkit.html
but have not tried it.

Agus

On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper <[hidden email]> wrote:
> Hi,
>
> I think you are asking for connected component labeling. raster::clump() expects your features to be like islands surrounded by NA or 0 background. There might be more Raster* centric ways to do it, but you could use plain old vanilla image processing using the very good imager package ( https://cran.r-project.org/web/packages/imager/ <https://cran.r-project.org/web/packages/imager/> ).  The steps below show in color your actual values with the connected component IDs labeled on each cell.  I used the 8-connected labeling but 4-connected is also available.
>
> library(raster)
> library(imager)
>
> set.seed(10)
> nc <- 8
> nr <- 8
> MIN <- 1
> MAX <- 4
> LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")
>
> r <- raster(ncols=nc, nrows=nr)
> r[] <- round(runif(ncell(r),MIN, MAX),digits=0)
> cc <- r
>
> img <- imager::as.cimg(as.matrix(r))
> labeled <- imager::label(img, high_connectivity = TRUE)
> cc[] <- as.matrix(labeled)
>
> plot(r, col = LUT)
> txt <- cc[]
> xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))
> text(xy[,1], xy[,2], txt )
>
>
> Cheers,
> Ben
>> On Jan 26, 2017, at 4:48 PM, Nelly Reduan <[hidden email]> wrote:
>>
>> Hello,
>> I would like to assign a unique ID number to each patch (a patch is composed of a set of adjacent cells) as shown in this figure:
>>  <2328.png>
>>
>> I tested the clump function (package raster) by applying an eight adjacent cells rule but this function mixes all cell values into single patch.
>>
>> Here is a code example to create a raster. The raster contains values ranged from 1 to 9.
>> r <- raster(ncols=12, nrows=12)
>> r[] <- round(runif(ncell(r),1,9),digits=0)
>> plot(r)
>>
>> Is there a way to assign a unique ID number to grouping cells that form a patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?
>>
>> Thanks a lot for your time
>> Nell
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email] <mailto:[hidden email]>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

raster brick from netcdf - specify range

Sat, 02/25/2017 - 13:03
Want to read part of a netcdf file into a raster brick. The file (in
ncdf4 order) contains a variable
pre[lon, lat, time]
and those dimensions are: 720, 360, 1356

The following reads in all time layers:
rb <- brick( filename, varname=varname, nl=2 )

> rb
class       : RasterBrick
dimensions  : 360, 720, 259200, 1356  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
z-value     : 12, 1367 (min, max)
varname     : pre

Likewise, I've had no success limiting the read using xmn, xmx, ymn,
ymx, or template. The help page for brick suggests that all these should
work, or at least, it does not seem say that they don't work for netcdf
files.

My workaround requires directly reading the netcdf file, then making a
brick from that array, but when reading finer resolution data, the
memory use and speed could be an issue due to creating 2 objects instead
of 1.

raster version = 2.5-8

I'd appreciate any suggestions about this.

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: split SLDF with attributes

Fri, 02/24/2017 - 15:13
Marta,
Since you have intersecting lines in your SLDF, another approach might
be to use a loop (!) through your SpatialLines objects and carry the
attributes, something like:

library(raster)

sldf_in <- shapefile("mergedDF_1_DF_short")

# Init a list to keep all segments
sldf_list <- vector("list", nrow(sldf_in))
names(sldf_list) <- row.names(sldf_in)

for (i in names(sldf_list)) {
   sl <- SegmentSpatialLines(as(sldf_in[i,], "SpatialLines"),
length=10000, merge.last=TRUE)
   df <- sldf_in[i,]@data
   df <- df[rep(seq_len(nrow(df)), length(sl)), ]
   df$ID <- i
   sldf_list[[i]] <- SpatialLinesDataFrame(sl, df, match.ID=FALSE)
   # Make arbitrary unique IDs (for combining in next step)
   sldf_list[[i]] <- spChFIDs(sldf_list[[i]], paste(i,
row.names(sldf_list[[i]]), sep="-"))
}

# Combine into a unique SLDF
sldf_out <- do.call(rbind, sldf_list)


--Mel.



On 2/24/2017 10:31 AM, marta azores wrote:
> Thanks for your answer Mel,
> Option1) Createsegments: this script cut the SLDF...but as autoridades
> said um the website...The results are spatiallines, not SLDF. Then,
> they don't jeep the attributes.
> Option2) rgeos:over I know this tool, but O have several lines
> crossing among them. I am afraid, that can be a bit messy.
> I doesn't know the gTouches...I check it after your suggestion..
> But I think the rgeos:Over os the best option for me.
> Thanks for your suggestion
> Marta
>
>
>
> No dia sexta-feira, 24 de fevereiro de 2017, Bacou, Melanie
> <[hidden email] <mailto:[hidden email]>> escreveu:
>
>     Marta,
>     Option 1) using CreateSegments() seems to do what you want. Why
>     don't you simply use rgeos::over() or rgeos::gTouches() between
>     your old and new SLDFs to match the index positions of your old
>     attributes to your new segments?
>
>     --Mel.
>
>
>     On 2/23/2017 2:27 PM, marta azores wrote:
>
>         Hi all,
>
>         I am working with a large datasets of boat positions, and
>         based on this
>         points I created a SpatialLinesDataFrame. Now, I would like to
>         split/divide SpatialLinesDataframe
>         into n segments of 10km and keep the attributes. So far, I
>         have tried
>         multiple things:
>
>         1.      The function from this website:
>         http://rstudio-pubs-static.s3.
>         amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html
>         <http://amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html>.
>         The problem with
>         it is that the function cuts the lines in segments of 10km but
>         it is not
>         possible to get the attributes.
>
>         2.      I have also tried the function from this website:
>         http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline-
>         <http://stackoverflow.com/questions/38700246/how-do-i-split-divide-polyline->
>         shapefiles-into-equally-length-smaller-segments. The problem
>         is that do not
>         cut the segments in 10km, and it seems to cut the segments
>         with random
>         lengths.
>
>         3.      On the follwoing website: http://r-sig-geo.2731867.n2.
>         nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#
>         <http://nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#>
>         a7583629 propose to divide lines by segments but do not
>         specify how it is
>         possible to do in based on the length of the segment.
>
>         Is there any function to do this in R? Or does somebody know
>         the solution
>         for that?
>
>         My data looks like this:
>
>         #Script:
>
>         path_data <- "C:/R/"
>
>         mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")
>
>
>         https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing
>         <https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklpd3ZnRzA?usp=sharing>
>
>                 [[alternative HTML version deleted]]
>
>         _______________________________________________
>         R-sig-Geo mailing list
>         [hidden email]
>         https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>         <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: split SLDF with attributes

Fri, 02/24/2017 - 09:31
Thanks for your answer Mel,
Option1) Createsegments: this script cut the SLDF...but as autoridades said
um the website...The results are spatiallines, not SLDF. Then, they don't
jeep the attributes.
Option2) rgeos:over I know this tool, but O have several lines crossing
among them. I am afraid, that can be a bit messy.
I doesn't know the gTouches...I check it after your suggestion..
But I think the rgeos:Over os the best option for me.
Thanks for your suggestion
Marta



No dia sexta-feira, 24 de fevereiro de 2017, Bacou, Melanie <[hidden email]>
escreveu:

> Marta,
> Option 1) using CreateSegments() seems to do what you want. Why don't you
> simply use rgeos::over() or rgeos::gTouches() between your old and new
> SLDFs to match the index positions of your old attributes to your new
> segments?
>
> --Mel.
>
>
> On 2/23/2017 2:27 PM, marta azores wrote:
>
>> Hi all,
>>
>> I am working with a large datasets of boat positions, and based on this
>> points I created a SpatialLinesDataFrame. Now, I would like to
>> split/divide SpatialLinesDataframe
>> into n segments of 10km and keep the attributes. So far, I have tried
>> multiple things:
>>
>> 1.      The function from this website: http://rstudio-pubs-static.s3.
>> amazonaws.com/10685_1f7266d60db7432486517a111c76ac8b.html. The problem
>> with
>> it is that the function cuts the lines in segments of 10km but it is not
>> possible to get the attributes.
>>
>> 2.      I have also tried the function from this website:
>> http://stackoverflow.com/questions/38700246/how-do-i-split-
>> divide-polyline-
>> shapefiles-into-equally-length-smaller-segments. The problem is that do
>> not
>> cut the segments in 10km, and it seems to cut the segments with random
>> lengths.
>>
>> 3.      On the follwoing website: http://r-sig-geo.2731867.n2.
>> nabble.com/split-divide-SpatialLines-sp-into-n-segments-td7583234.html#
>> a7583629 propose to divide lines by segments but do not specify how it is
>> possible to do in based on the length of the segment.
>>
>> Is there any function to do this in R? Or does somebody know the solution
>> for that?
>>
>> My data looks like this:
>>
>> #Script:
>>
>> path_data <- "C:/R/"
>>
>> mergedDF_1_DF_short<- shapefile(path_data,"mergedDF_1_DF_short")
>>
>>
>> https://drive.google.com/drive/folders/0B7IbvWhE5JNPT0oxUklp
>> d3ZnRzA?usp=sharing
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Pages