This is an

**for**__archive__**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

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

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

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

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?

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

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?

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

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

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

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

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

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

*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

*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

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

> 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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