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:*4 min 11 sec ago

### Making an hexagonal grid using spsample

Dear list members,

I am trying to make an spatial hexagonal grid in R using the spsample

function. I would like hexagons of 1 ha = 10000 sq meters.

First I used a cell size of 62.04 m as the side of the hexagon.

data(meuse.grid)

gridded(meuse.grid) = ~x + y

proj4string(meuse.grid) <- CRS("+init=epsg:28992")

xx <- spsample(meuse.grid, type = "hexagonal", cellsize = 62.04)

xx <- HexPoints2SpatialPolygons(xx)

rgeos::gArea(xx, byid=TRUE)

It gave me an area per hexagon = 3333.299

Second I used a cell size of 62.04 * 2 = 124.08 m as the long diagonal of

the hexagon

data(meuse.grid)

gridded(meuse.grid) = ~x + y

proj4string(meuse.grid) <- CRS("+init=epsg:28992")

xx <- spsample(meuse.grid, type = "hexagonal", cellsize = 62.04*2)

xx <- HexPoints2SpatialPolygons(xx)

rgeos::gArea(xx, byid=TRUE)

It gave me an area per hexagon = 13333.19

Both results are wrong, because I was expecting close to 10000.

How can I make a spatial grid of 1 ha (or other size) in R?

Manuel

--

*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

I am trying to make an spatial hexagonal grid in R using the spsample

function. I would like hexagons of 1 ha = 10000 sq meters.

First I used a cell size of 62.04 m as the side of the hexagon.

data(meuse.grid)

gridded(meuse.grid) = ~x + y

proj4string(meuse.grid) <- CRS("+init=epsg:28992")

xx <- spsample(meuse.grid, type = "hexagonal", cellsize = 62.04)

xx <- HexPoints2SpatialPolygons(xx)

rgeos::gArea(xx, byid=TRUE)

It gave me an area per hexagon = 3333.299

Second I used a cell size of 62.04 * 2 = 124.08 m as the long diagonal of

the hexagon

data(meuse.grid)

gridded(meuse.grid) = ~x + y

proj4string(meuse.grid) <- CRS("+init=epsg:28992")

xx <- spsample(meuse.grid, type = "hexagonal", cellsize = 62.04*2)

xx <- HexPoints2SpatialPolygons(xx)

rgeos::gArea(xx, byid=TRUE)

It gave me an area per hexagon = 13333.19

Both results are wrong, because I was expecting close to 10000.

How can I make a spatial grid of 1 ha (or other size) in R?

Manuel

--

*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: Create Nodes at lines' segments extremities.

(after receiving the data, off-list):

library(sp)

library(rgdal)

x = readOGR("TRONCON_TOPO_Bresle.shp")

lapply(x@lines, function(y) lapply(y@Lines, function(z) apply(z@coords,

2, range)))

gives you the bounding box of each individual Line element, as a nested

list (similar to how x@lines is nested).

On 02/03/17 17:53, Tristan Bourgeois wrote:

> Hi everybody,

>

> I can't find any infomation about the gis manipulation I want to process on

> R. Even if there are similar posts on r-sig-geo archives...

>

>

> My issue is the following one :

>

>> *myline<-readOGR(dsn = "./shapes/shp_tr_topo",layer =

> "TRONCON_TOPO_Bresle" ) is a spatial polyline which is composed of :*

>

>> *length(line@lines)*

> [1] 585 => 585 segments

>

>

> It seems that each of those 585 lines are acually a serie of points :

>

> *> line@lines[[1]]@Lines*

> [[1]]

> An object of class "Line"

> Slot "coords":

> [,1] [,2]

> [1,] 609046.7 6956543

> [2,] 609057.3 6956543

> [3,] 609077.6 6956546

> [4,] 609103.2 6956553

> [5,] 609132.5 6956561

> [6,] 609146.0 6956568

> [7,] 609157.4 6956574

> [8,] 609167.0 6956577

> [9,] 609176.9 6956581

> [10,] 609184.4 6956585

> [11,] 609191.9 6956588

> [12,] 609199.3 6956591

> [13,] 609205.7 6956594

> [14,] 609212.6 6956599

> [15,] 609216.6 6956604

> [16,] 609222.6 6956609

> [17,] 609229.0 6956613

> [18,] 609237.9 6956619

> [19,] 609244.4 6956623

> [20,] 609253.5 6956626

> [21,] 609265.2 6956631

> [22,] 609276.5 6956636

> [23,] 609287.4 6956639

> [24,] 609296.9 6956644

> [25,] 609305.1 6956647

> [26,] 609308.6 6956653

>

> What I want to obtain are the coordinates of each segments extremities. So

> for a segment I should have 2 points with X,Y values.

> I tried to apply the rgeos Gnode function but nothing appears as changed.

> Maybe there's something to do with the min and max X,Y of each segments ?

>

>

> Does somebody has already done /wrote a similar script to process that ??

>

>

> Thank you in advance for the helpfull advices.

>

> Cheers.

>

> Tristan Bourgeois

> Hydromorphologist Engineer

> Seine Normandie Water Agency

> Rouen, France

>

> [[alternative HTML version deleted]]

>

> _______________________________________________

> R-sig-Geo mailing list

> [hidden email]

> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

> --

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: Suggestions for irregular grid netCDF data?

Hi again,

Just to follow-up - got this working without too much issue... as of now, I

have the data values associated w/ the center lat/long of each grid cell.

There are also lat/long vertices for polygons associated with each grid

cell, but that'll take more working through it to get it together. However,

for now this largely suits my needs, and for now I can approximate the grid

cells using voronoi polygons for the grid cell centroids.

Here's some working code [I stop at creating dataframes, but going the next

steps for spatial objects is pretty easy - either sp or simple features

(still learning the latter)].

Thanks again!

Mike

#Sample Code

#Gets netCDF data for lat/long points

library(ncdf4)

#Connect to netCDF

nyhops.nc <- nc_open("

http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc

")

##To get 'simple', 3D variale (e.g., depth)

#Load depth, lat, and long

depth.nc <- ncvar_get(nyhops.nc, 'depth')

lat.nc <- ncvar_get(nyhops.nc, 'lat')

long.nc <- ncvar_get(nyhops.nc, 'lon')

#Create dataframe

depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),

depth=as.vector(depth.nc))

#Remove NA values

depth.df <- depth.df[!is.na(depth.df$depth),]

#Can then convert to spatial object, export as table, etc.

##Example getting data involving a vertical horizon dimension and a date

dimension

#Fetch the temperature data; start is the starting position for x, y,

sigma, date; count is how many

# here, all rows & columns, and one horizon (horizon 3) and date (394)

temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(1,1,3,394),

count=c(147,452,1,1))

#Can Find out how many dimensions and entries:

nyhops.nc$var[['temp']]$varsize

#Output: [1] 147 452 11 394

#147 rows; 452 columns; 11 horizons; 394 dates/times

#Create dataframe w/ lat, long, temp

temp.df <- data.frame(cbind(long=as.vector(long.nc), lat=as.vector(lat.nc),

temp=as.vector(temp.nc)))

#remove NA rows for temp

temp.df <- temp.df[!is.na(temp.df$temp),] #Can conver to spatial object,

export as table, etc.

##Can plot, just as a non-spatial object for now.

#assign color ramp

rbPal <- colorRampPalette(c('blue', 'red'))

#Adds a column of color values based on the temp values

temp.df$Col <- rbPal(10)[as.numeric(cut(temp.df$temp,breaks = 10))]

plot(temp.df$long, temp.df$lat, col=temp.df$Col)

##Next step - use lon_vertices and lat_vertices to get polygons associated

with each point

On Thu, Mar 2, 2017 at 8:58 AM, Michael Treglia <[hidden email]> wrote:

> Thanks Mike!

>

> This is helpful, so far at least to see/think about different ways to work

> with the data, as I wade through this structure and understand how

> different R packages handle it.

>

> I'll give things a go, and follow up, and will keep an eye out for any

> other suggestions. And if you catch anything I'm missing if you go through

> the example code, I'm all ears. This isn't pressing on my end, but working

> through this dataset (and similarly structured ones) will ultimately be

> really useful.

>

> Thanks again,

> Mike T

>

> Please pardon any typos, this message was sent from a mobile device.

>

> On Mar 1, 2017 11:43 PM, "Michael Sumner" <[hidden email]> wrote:

>

>> 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/Hindc

>>> ast/all_nb_mon_81.nc,

>>> and here’s the web-interface for the data:

>>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindc

>>> ast/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/Hindc

>>> ast/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/Hindc

>>> ast/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

Just to follow-up - got this working without too much issue... as of now, I

have the data values associated w/ the center lat/long of each grid cell.

There are also lat/long vertices for polygons associated with each grid

cell, but that'll take more working through it to get it together. However,

for now this largely suits my needs, and for now I can approximate the grid

cells using voronoi polygons for the grid cell centroids.

Here's some working code [I stop at creating dataframes, but going the next

steps for spatial objects is pretty easy - either sp or simple features

(still learning the latter)].

Thanks again!

Mike

#Sample Code

#Gets netCDF data for lat/long points

library(ncdf4)

#Connect to netCDF

nyhops.nc <- nc_open("

http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc

")

##To get 'simple', 3D variale (e.g., depth)

#Load depth, lat, and long

depth.nc <- ncvar_get(nyhops.nc, 'depth')

lat.nc <- ncvar_get(nyhops.nc, 'lat')

long.nc <- ncvar_get(nyhops.nc, 'lon')

#Create dataframe

depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),

depth=as.vector(depth.nc))

#Remove NA values

depth.df <- depth.df[!is.na(depth.df$depth),]

#Can then convert to spatial object, export as table, etc.

##Example getting data involving a vertical horizon dimension and a date

dimension

#Fetch the temperature data; start is the starting position for x, y,

sigma, date; count is how many

# here, all rows & columns, and one horizon (horizon 3) and date (394)

temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(1,1,3,394),

count=c(147,452,1,1))

#Can Find out how many dimensions and entries:

nyhops.nc$var[['temp']]$varsize

#Output: [1] 147 452 11 394

#147 rows; 452 columns; 11 horizons; 394 dates/times

#Create dataframe w/ lat, long, temp

temp.df <- data.frame(cbind(long=as.vector(long.nc), lat=as.vector(lat.nc),

temp=as.vector(temp.nc)))

#remove NA rows for temp

temp.df <- temp.df[!is.na(temp.df$temp),] #Can conver to spatial object,

export as table, etc.

##Can plot, just as a non-spatial object for now.

#assign color ramp

rbPal <- colorRampPalette(c('blue', 'red'))

#Adds a column of color values based on the temp values

temp.df$Col <- rbPal(10)[as.numeric(cut(temp.df$temp,breaks = 10))]

plot(temp.df$long, temp.df$lat, col=temp.df$Col)

##Next step - use lon_vertices and lat_vertices to get polygons associated

with each point

On Thu, Mar 2, 2017 at 8:58 AM, Michael Treglia <[hidden email]> wrote:

> Thanks Mike!

>

> This is helpful, so far at least to see/think about different ways to work

> with the data, as I wade through this structure and understand how

> different R packages handle it.

>

> I'll give things a go, and follow up, and will keep an eye out for any

> other suggestions. And if you catch anything I'm missing if you go through

> the example code, I'm all ears. This isn't pressing on my end, but working

> through this dataset (and similarly structured ones) will ultimately be

> really useful.

>

> Thanks again,

> Mike T

>

> Please pardon any typos, this message was sent from a mobile device.

>

> On Mar 1, 2017 11:43 PM, "Michael Sumner" <[hidden email]> wrote:

>

>> 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/Hindc

>>> ast/all_nb_mon_81.nc,

>>> and here’s the web-interface for the data:

>>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindc

>>> ast/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/Hindc

>>> ast/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/Hindc

>>> ast/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

### Create Nodes at lines' segments extremities.

Hi everybody,

I can't find any infomation about the gis manipulation I want to process on

R. Even if there are similar posts on r-sig-geo archives...

My issue is the following one :

>*myline<-readOGR(dsn = "./shapes/shp_tr_topo",layer =

"TRONCON_TOPO_Bresle" ) is a spatial polyline which is composed of :*

> *length(line@lines)*

[1] 585 => 585 segments

It seems that each of those 585 lines are acually a serie of points :

*> line@lines[[1]]@Lines*

[[1]]

An object of class "Line"

Slot "coords":

[,1] [,2]

[1,] 609046.7 6956543

[2,] 609057.3 6956543

[3,] 609077.6 6956546

[4,] 609103.2 6956553

[5,] 609132.5 6956561

[6,] 609146.0 6956568

[7,] 609157.4 6956574

[8,] 609167.0 6956577

[9,] 609176.9 6956581

[10,] 609184.4 6956585

[11,] 609191.9 6956588

[12,] 609199.3 6956591

[13,] 609205.7 6956594

[14,] 609212.6 6956599

[15,] 609216.6 6956604

[16,] 609222.6 6956609

[17,] 609229.0 6956613

[18,] 609237.9 6956619

[19,] 609244.4 6956623

[20,] 609253.5 6956626

[21,] 609265.2 6956631

[22,] 609276.5 6956636

[23,] 609287.4 6956639

[24,] 609296.9 6956644

[25,] 609305.1 6956647

[26,] 609308.6 6956653

What I want to obtain are the coordinates of each segments extremities. So

for a segment I should have 2 points with X,Y values.

I tried to apply the rgeos Gnode function but nothing appears as changed.

Maybe there's something to do with the min and max X,Y of each segments ?

Does somebody has already done /wrote a similar script to process that ??

Thank you in advance for the helpfull advices.

Cheers.

Tristan Bourgeois

Hydromorphologist Engineer

Seine Normandie Water Agency

Rouen, France

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

I can't find any infomation about the gis manipulation I want to process on

R. Even if there are similar posts on r-sig-geo archives...

My issue is the following one :

>*myline<-readOGR(dsn = "./shapes/shp_tr_topo",layer =

"TRONCON_TOPO_Bresle" ) is a spatial polyline which is composed of :*

> *length(line@lines)*

[1] 585 => 585 segments

It seems that each of those 585 lines are acually a serie of points :

*> line@lines[[1]]@Lines*

[[1]]

An object of class "Line"

Slot "coords":

[,1] [,2]

[1,] 609046.7 6956543

[2,] 609057.3 6956543

[3,] 609077.6 6956546

[4,] 609103.2 6956553

[5,] 609132.5 6956561

[6,] 609146.0 6956568

[7,] 609157.4 6956574

[8,] 609167.0 6956577

[9,] 609176.9 6956581

[10,] 609184.4 6956585

[11,] 609191.9 6956588

[12,] 609199.3 6956591

[13,] 609205.7 6956594

[14,] 609212.6 6956599

[15,] 609216.6 6956604

[16,] 609222.6 6956609

[17,] 609229.0 6956613

[18,] 609237.9 6956619

[19,] 609244.4 6956623

[20,] 609253.5 6956626

[21,] 609265.2 6956631

[22,] 609276.5 6956636

[23,] 609287.4 6956639

[24,] 609296.9 6956644

[25,] 609305.1 6956647

[26,] 609308.6 6956653

What I want to obtain are the coordinates of each segments extremities. So

for a segment I should have 2 points with X,Y values.

I tried to apply the rgeos Gnode function but nothing appears as changed.

Maybe there's something to do with the min and max X,Y of each segments ?

Does somebody has already done /wrote a similar script to process that ??

Thank you in advance for the helpfull advices.

Cheers.

Tristan Bourgeois

Hydromorphologist Engineer

Seine Normandie Water Agency

Rouen, France

[[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?

Thanks Mike!

This is helpful, so far at least to see/think about different ways to work

with the data, as I wade through this structure and understand how

different R packages handle it.

I'll give things a go, and follow up, and will keep an eye out for any

other suggestions. And if you catch anything I'm missing if you go through

the example code, I'm all ears. This isn't pressing on my end, but working

through this dataset (and similarly structured ones) will ultimately be

really useful.

Thanks again,

Mike T

Please pardon any typos, this message was sent from a mobile device.

On Mar 1, 2017 11:43 PM, "Michael Sumner" <[hidden email]> wrote:

> 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

This is helpful, so far at least to see/think about different ways to work

with the data, as I wade through this structure and understand how

different R packages handle it.

I'll give things a go, and follow up, and will keep an eye out for any

other suggestions. And if you catch anything I'm missing if you go through

the example code, I'm all ears. This isn't pressing on my end, but working

through this dataset (and similarly structured ones) will ultimately be

really useful.

Thanks again,

Mike T

Please pardon any typos, this message was sent from a mobile device.

On Mar 1, 2017 11:43 PM, "Michael Sumner" <[hidden email]> wrote:

> 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

### 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