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

Making an hexagonal grid using spsample

Fri, 03/03/2017 - 13:03
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

Re: Create Nodes at lines' segments extremities.

Fri, 03/03/2017 - 08:09

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

Thu, 03/02/2017 - 22:11
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

Create Nodes at lines' segments extremities.

Thu, 03/02/2017 - 10:53
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

Re: Suggestions for irregular grid netCDF data?

Thu, 03/02/2017 - 07:58
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

ApproxNA producing different and incorrect results in stacks of large rasters

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



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



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



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



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



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



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



Best regards, and thanks for any suggestions



Jodi



##### start  Reproducible example



library(raster)

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

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

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

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

r2 <- setValues(r, NA)

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

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

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

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



# insert some NA values

r1[100:600,] <- NA

r3[,1500:1950] <- NA

r5[,400:800] <- NA

r6[2500:2900,] <- NA



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

r4[300:500,] <- 750

r6[300:400,] <- 20

r6[400:500,] <- 100

# make a stack

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



# see what the pre-interpolation stack looks like

plot(s)



# interpolate, this takes a while

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



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

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

plot(x2)



## End reproducible example







--please do not edit the information below--

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

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

Windows 7 x64 (build 7601) Service Pack 1

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

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







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


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

        [[alternative HTML version deleted]]

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

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

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

Manuel

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

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


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

        [[alternative HTML version deleted]]

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

Re: Suggestions for irregular grid netCDF data?

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

my general advice is

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

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

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

http://rpubs.com/cyclemumner/roms0

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

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

Cheers, Mike.


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

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

        [[alternative HTML version deleted]]

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

Suggestions for irregular grid netCDF data?

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

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

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

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

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

Thanks!
Mike


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

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

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

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

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



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

        [[alternative HTML version deleted]]

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

SpatialDataFrame List Union

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

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

My code :

library(rgdal)

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

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


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

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

Does anybody has a pist ?

Thanks !

Cheers

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

        [[alternative HTML version deleted]]

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

EGU: R’s deliberate role in Earth sciences

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

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


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

Re: split SLDF with attributes

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

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


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


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

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

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

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

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

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

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

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




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

Thanks in advance,

Marta

>
>>
>

        [[alternative HTML version deleted]]

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

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

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

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

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

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

cheers,

Rolf Turner

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

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

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

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

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


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


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

        [[alternative HTML version deleted]]

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

Re: rasterTopolygons: singleparts often needed

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

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

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

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

it works for me:

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

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

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


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

        [[alternative HTML version deleted]]

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

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

Sun, 02/26/2017 - 11:46


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

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

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

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

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

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

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


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

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

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

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

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

plot(nc["SID79"])



Best,

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

        [[alternative HTML version deleted]]

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

Re: rasterTopolygons: singleparts often needed

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

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

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

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

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

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

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

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

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

Pages