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: 8 min 30 sec ago

BIG DATABASE

Thu, 05/24/2018 - 06:29
Hello,

Is it possible with R to work on a big database without loading it? If yes,
how can I do do it?

Thanks.

--
Yaya BAMBA

Elève Ingénieur Statisticien Economiste (ISE)

Ecole Nationale Supérieure de Statistique et d'Economie Appliquée (ENSEA),
Abidjan (Côte d'Ivoire)

Tél: +225 87 89 76 89

        [[alternative HTML version deleted]]

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

Re: how to plot different rows of a SpatialPolygonsDataFrame in trellis panels

Wed, 05/23/2018 - 23:28
Dear Scott,

use
spplot(tmp, zcol="z")
instead of
spplot(tmp, zcol=z)
to display the column "z".

Use
spplot(tmp[tmp$id == 1, ], zcol = "z")
to display only those grids that has id 1.

And to display two panels:
tmp$z1 <- tmp$z2 <- tmp$z
tmp$z1[tmp$id == 1] <-  NA
tmp$z2[tmp$id == 2] <-  NA
spplot(tmp, zcol=c("z1", "z2"))

HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences


2018.05.23. 21:50 keltezéssel, Waichler, Scott R írta:
> Roger,
>
>> Without having tried your code, how similar is this to these packages:
>>
>> https://cran.r-project.org/package=micromap
>> https://cran.r-project.org/package=micromapST
> Thanks for responding.  I looked at the micromap vignettes but it seems their focus is being able to provide small maps in one column in figures where the main emphasis is on statistics presented in other columns.  In my problem, all I want to show are maps.  I haven't been able to do this with a higher-level function like spplot or stplot.  I can imagine doing it with a custom panel function, where lpolygons() is used with input that depends on the panel.
>
> In case it helps, I should clarify that the "different polygons" I referred to are polygons with different spatial coordinates; the attribute values are a small number of discrete values.  My polygons are representations of contaminant plumes in groundwater.  For each year, I have read in a shapefile containing polygons of the contaminant concentration at 2 to 4 levels.  I have combined all of these data into a SpatialPolygonsDataFrame, where each row represents a given year and concentration level, and contains one or more polygons.  The coordinates are different in each row, but years and concentration level are common across rows.  This seems to match the nature of a "long table" format in spacetime, but I haven't been able to work it out in that package, hence I'm stepping back and trying to get the basics working with spplot.
>
> Thanks,
> Scott
>
>>> Hello,
>>>
>>> I have a SpatialPolygonsDataFrame.  I would like to do a trellis plot on one
>> of the attributes, so that in the panel for a given attribute value, only those
>> polygons with that value are plotted.  So, each panel has different polygons
>> plotted in it.  I can't figure out how to do this.  In the toy example below, I
>> would like to create a trellis plot with one panel showing the polygons with id
>> = 1, and another panel showing the polygons with id = 2.
>>> My goal beyond this toy problem is to do the same thing with stplot, where
>> panels correspond to times and each time has a different set of polygons
>> plotted.  Will that be possible?  In all the examples I can find of using stplot
>> for a space-time grid with the spatial objects being polygons, the polygons
>> are the same across time.
>>> # based on example in help("SpatialPolygonsDataFrame-class")
>>> Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
>>> Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
>>> Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
>>> Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
>>> Srs1 = Polygons(list(Sr1), "s1")
>>> Srs2 = Polygons(list(Sr2), "s2")
>>> Srs3 = Polygons(list(Sr3, Sr4), "s3/4") SpP =
>>> SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3) grd <- GridTopology(c(1,1),
>>> c(1,1), c(10,10)) polys <- as(grd, "SpatialPolygons") centroids <-
>>> coordinates(polys) x <- centroids[,1] y <- centroids[,2] z <- 1.4 +
>>> 0.1*x + 0.2*y + 0.002*x*x id = factor(sample(c(1,2),
>>> size=length(polys), replace=T)) tmp <- SpatialPolygonsDataFrame(polys,
>>>       data=data.frame(x=x, y=y, z=z, id=id,
>>> row.names=row.names(polys)))
>>> plot(tmp)  # plots all the square polygons (n = 10*10)
>>> spplot(tmp)  # plots values of x, y, z, id in separate panels, each
>>> with 100 polys spplot(tmp, zcol=z)  # error message about duplication
>>> of factor level spplot(tmp ~ id, zcol=z, data=tmp)  # won't take
>>> formula
>>>
>>> Thank you,
>>> ScottWaichler
>>> Pacific Northwest National Laboratory
>>> scott.waichler _at_ pnnl.gov
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics, Helleveien 30, N-
>> 5045 Bergen, Norway.
>> voice: +47 55 95 93 55; e-mail: [hidden email]
>> http://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
> _______________________________________________
> 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 plot different rows of a SpatialPolygonsDataFrame in trellis panels

Wed, 05/23/2018 - 14:50
Roger,

> Without having tried your code, how similar is this to these packages:
>
> https://cran.r-project.org/package=micromap
> https://cran.r-project.org/package=micromapST

Thanks for responding.  I looked at the micromap vignettes but it seems their focus is being able to provide small maps in one column in figures where the main emphasis is on statistics presented in other columns.  In my problem, all I want to show are maps.  I haven't been able to do this with a higher-level function like spplot or stplot.  I can imagine doing it with a custom panel function, where lpolygons() is used with input that depends on the panel.  

In case it helps, I should clarify that the "different polygons" I referred to are polygons with different spatial coordinates; the attribute values are a small number of discrete values.  My polygons are representations of contaminant plumes in groundwater.  For each year, I have read in a shapefile containing polygons of the contaminant concentration at 2 to 4 levels.  I have combined all of these data into a SpatialPolygonsDataFrame, where each row represents a given year and concentration level, and contains one or more polygons.  The coordinates are different in each row, but years and concentration level are common across rows.  This seems to match the nature of a "long table" format in spacetime, but I haven't been able to work it out in that package, hence I'm stepping back and trying to get the basics working with spplot.

Thanks,
Scott

> > Hello,
> >
> > I have a SpatialPolygonsDataFrame.  I would like to do a trellis plot on one
> of the attributes, so that in the panel for a given attribute value, only those
> polygons with that value are plotted.  So, each panel has different polygons
> plotted in it.  I can't figure out how to do this.  In the toy example below, I
> would like to create a trellis plot with one panel showing the polygons with id
> = 1, and another panel showing the polygons with id = 2.
> >
> > My goal beyond this toy problem is to do the same thing with stplot, where
> panels correspond to times and each time has a different set of polygons
> plotted.  Will that be possible?  In all the examples I can find of using stplot
> for a space-time grid with the spatial objects being polygons, the polygons
> are the same across time.
> >
> > # based on example in help("SpatialPolygonsDataFrame-class")
> > Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
> > Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
> > Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
> > Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
> > Srs1 = Polygons(list(Sr1), "s1")
> > Srs2 = Polygons(list(Sr2), "s2")
> > Srs3 = Polygons(list(Sr3, Sr4), "s3/4") SpP =
> > SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3) grd <- GridTopology(c(1,1),
> > c(1,1), c(10,10)) polys <- as(grd, "SpatialPolygons") centroids <-
> > coordinates(polys) x <- centroids[,1] y <- centroids[,2] z <- 1.4 +
> > 0.1*x + 0.2*y + 0.002*x*x id = factor(sample(c(1,2),
> > size=length(polys), replace=T)) tmp <- SpatialPolygonsDataFrame(polys,
> >      data=data.frame(x=x, y=y, z=z, id=id,
> > row.names=row.names(polys)))
> > plot(tmp)  # plots all the square polygons (n = 10*10)
> > spplot(tmp)  # plots values of x, y, z, id in separate panels, each
> > with 100 polys spplot(tmp, zcol=z)  # error message about duplication
> > of factor level spplot(tmp ~ id, zcol=z, data=tmp)  # won't take
> > formula
> >
> > Thank you,
> > ScottWaichler
> > Pacific Northwest National Laboratory
> > scott.waichler _at_ pnnl.gov
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics, Helleveien 30, N-
> 5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: how to plot different rows of a SpatialPolygonsDataFrame in trellis panels

Wed, 05/23/2018 - 13:53
Scott,

Without having tried your code, how similar is this to these packages:

https://cran.r-project.org/package=micromap
https://cran.r-project.org/package=micromapST

and their JSS papers:

https://www.jstatsoft.org/index.php/jss/article/view/v063i02
https://www.jstatsoft.org/index.php/jss/article/view/v063i03

Roger

On Wed, 23 May 2018, Waichler, Scott R wrote:

> Hello,
>
> I have a SpatialPolygonsDataFrame.  I would like to do a trellis plot on one of the attributes, so that in the panel for a given attribute value, only those polygons with that value are plotted.  So, each panel has different polygons plotted in it.  I can't figure out how to do this.  In the toy example below, I would like to create a trellis plot with one panel showing the polygons with id = 1, and another panel showing the polygons with id = 2.
>
> My goal beyond this toy problem is to do the same thing with stplot, where panels correspond to times and each time has a different set of polygons plotted.  Will that be possible?  In all the examples I can find of using stplot for a space-time grid with the spatial objects being polygons, the polygons are the same across time.
>
> # based on example in help("SpatialPolygonsDataFrame-class")
> Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
> Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
> Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
> Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
> Srs1 = Polygons(list(Sr1), "s1")
> Srs2 = Polygons(list(Sr2), "s2")
> Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
> SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
> grd <- GridTopology(c(1,1), c(1,1), c(10,10))
> polys <- as(grd, "SpatialPolygons")
> centroids <- coordinates(polys)
> x <- centroids[,1]
> y <- centroids[,2]
> z <- 1.4 + 0.1*x + 0.2*y + 0.002*x*x
> id = factor(sample(c(1,2), size=length(polys), replace=T))
> tmp <- SpatialPolygonsDataFrame(polys,
>      data=data.frame(x=x, y=y, z=z, id=id, row.names=row.names(polys)))
> plot(tmp)  # plots all the square polygons (n = 10*10)
> spplot(tmp)  # plots values of x, y, z, id in separate panels, each with 100 polys
> spplot(tmp, zcol=z)  # error message about duplication of factor level
> spplot(tmp ~ id, zcol=z, data=tmp)  # won't take formula
>
> Thank you,
> ScottWaichler
> Pacific Northwest National Laboratory
> scott.waichler _at_ pnnl.gov
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

how to plot different rows of a SpatialPolygonsDataFrame in trellis panels

Wed, 05/23/2018 - 13:39
Hello,

I have a SpatialPolygonsDataFrame.  I would like to do a trellis plot on one of the attributes, so that in the panel for a given attribute value, only those polygons with that value are plotted.  So, each panel has different polygons plotted in it.  I can't figure out how to do this.  In the toy example below, I would like to create a trellis plot with one panel showing the polygons with id = 1, and another panel showing the polygons with id = 2.

My goal beyond this toy problem is to do the same thing with stplot, where panels correspond to times and each time has a different set of polygons plotted.  Will that be possible?  In all the examples I can find of using stplot for a space-time grid with the spatial objects being polygons, the polygons are the same across time.

# based on example in help("SpatialPolygonsDataFrame-class")
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
grd <- GridTopology(c(1,1), c(1,1), c(10,10))
polys <- as(grd, "SpatialPolygons")
centroids <- coordinates(polys)
x <- centroids[,1]
y <- centroids[,2]
z <- 1.4 + 0.1*x + 0.2*y + 0.002*x*x
id = factor(sample(c(1,2), size=length(polys), replace=T))
tmp <- SpatialPolygonsDataFrame(polys,
      data=data.frame(x=x, y=y, z=z, id=id, row.names=row.names(polys)))
plot(tmp)  # plots all the square polygons (n = 10*10)
spplot(tmp)  # plots values of x, y, z, id in separate panels, each with 100 polys
spplot(tmp, zcol=z)  # error message about duplication of factor level
spplot(tmp ~ id, zcol=z, data=tmp)  # won't take formula

Thank you,
ScottWaichler
Pacific Northwest National Laboratory
scott.waichler _at_ pnnl.gov

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

Re: spatio temporal anisotropy: help with stAni argument

Tue, 05/22/2018 - 03:02
This question has also been posted to gis.stackexchange.com:

https://gis.stackexchange.com/questions/283607/spatio-temporal-anisotropy-in-r-stani-argument-in-vgmst-and-krigest/

and better answers than mine there are welcome...

Barry


On Tue, May 22, 2018 at 8:39 AM, Dr. Benedikt Gräler <[hidden email]>
wrote:

> Dear Bruno,
>
> stAni relates spatial and temporal distances to each other in the model
> for the data set under study. A value of 100 km/day reflects that two
> stations 100 km apart have about the same strength of dependence as the
> values one day apart at the same location (This is in general only true for
> at most one pair of spatial and temporal distances and a 'global'
> compromise is sought). The function gstat::estiStAni provides a few
> heuristics to estimate the spatio-temporal anisotropy.
>
> When a spatio-temporal vgm including a metric component is fitted, the
> stAni parameter is optimized under the objective to find the theoretical
> vgm surface that best fits the empirical one.
>
> In krigeST, an additional and approximate stAni parameter can be provided
> to ease the (pre-)selection of the most relevant observations in the
> space-time cube for the prediction when local kriging is carried out. Note
> that the most relevant (i.e. strongest correlated) observations need not to
> form a cube, but depend on the actual vgm-model in use (see also the
> vignette "Spatio-Temporal Interpolation using gstat" at [1]).
>
> The unit of stAni depends on the spatial and temporal units of your data
> set. Check the units of the empirical variogram to find out which units
> R/gstat uses for your data set.
>
> Miss-specifying the stAni parameter will put too large or too small
> weights on the observations with a temporal difference resulting in bad
> predictions (this also holds vise-versa for spatial distances).
>
> HTH,
>
>  Ben
>
> [1] https://cran.r-project.org/web/packages/gstat/vignettes/spat
> io-temporal-kriging.pdf
>
>
> On 21.05.2018 15:16, Bruno Sesti wrote:
>
>> Hi,
>> I have some problems in understanding the usage of argument 'stAni' of
>> vgmST and krigeST functions, for.the setting of spatio temporal
>> anisotropy.
>> How have I to interpret this argument, that is the integer value that I
>> saw
>> I have to assign to it, what does it means?
>> How can this argument influence the result of prediction?
>> Are units of spatial and temporal ranges and lags related to this
>> argument?
>> How can I study the spatio temporal anisotropy of my data?
>> Kind regards.
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
> --
> Dr. Benedikt Gräler
> 52°North Initiative for Geospatial Open Source Software GmbH
> Martin-Luther-King-Weg 24
> 48155 Muenster, Germany
>
> E-Mail: [hidden email]
> Fon: +49-(0)-251/396371-39
> Fax: +49-(0)-251/396371-11
>
> http://52north.org/
> Twitter: @FiveTwoN
>
> General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
> Local Court Muenster HRB 10849
>
> _______________________________________________
> 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: spatio temporal anisotropy: help with stAni argument

Tue, 05/22/2018 - 02:39
Dear Bruno,

stAni relates spatial and temporal distances to each other in the model
for the data set under study. A value of 100 km/day reflects that two
stations 100 km apart have about the same strength of dependence as the
values one day apart at the same location (This is in general only true
for at most one pair of spatial and temporal distances and a 'global'
compromise is sought). The function gstat::estiStAni provides a few
heuristics to estimate the spatio-temporal anisotropy.

When a spatio-temporal vgm including a metric component is fitted, the
stAni parameter is optimized under the objective to find the theoretical
vgm surface that best fits the empirical one.

In krigeST, an additional and approximate stAni parameter can be
provided to ease the (pre-)selection of the most relevant observations
in the space-time cube for the prediction when local kriging is carried
out. Note that the most relevant (i.e. strongest correlated)
observations need not to form a cube, but depend on the actual vgm-model
in use (see also the vignette "Spatio-Temporal Interpolation using
gstat" at [1]).

The unit of stAni depends on the spatial and temporal units of your data
set. Check the units of the empirical variogram to find out which units
R/gstat uses for your data set.

Miss-specifying the stAni parameter will put too large or too small
weights on the observations with a temporal difference resulting in bad
predictions (this also holds vise-versa for spatial distances).

HTH,

  Ben

[1]
https://cran.r-project.org/web/packages/gstat/vignettes/spatio-temporal-kriging.pdf


On 21.05.2018 15:16, Bruno Sesti wrote:
> Hi,
> I have some problems in understanding the usage of argument 'stAni' of
> vgmST and krigeST functions, for.the setting of spatio temporal anisotropy.
> How have I to interpret this argument, that is the integer value that I saw
> I have to assign to it, what does it means?
> How can this argument influence the result of prediction?
> Are units of spatial and temporal ranges and lags related to this argument?
> How can I study the spatio temporal anisotropy of my data?
> Kind regards.
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Dr. Benedikt Gräler
52°North Initiative for Geospatial Open Source Software GmbH
Martin-Luther-King-Weg 24
48155 Muenster, Germany

E-Mail: [hidden email]
Fon: +49-(0)-251/396371-39
Fax: +49-(0)-251/396371-11

http://52north.org/
Twitter: @FiveTwoN

General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
Local Court Muenster HRB 10849

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

stConstruct with multiple SpatialPolygonsDataFrame objects

Mon, 05/21/2018 - 19:40
Hi, Can anyone provide an example of using stConstruct() with multiple SpatialPolygonsDataFrame objects?  I have created the latter by reading in multiple shapefiles, one shapefile per year.  Each shapefile contains polygons at several levels of an attribute.  I want to plot these polygons for levels of say 1, 5, 10  using stplot, where the panels are years.  I am having trouble understanding how to do this type of task where the spatial coordinates vary across time.

Thank you,
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler _at_ pnnl.gov

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

Running spacetime eof on large raster stack

Mon, 05/21/2018 - 18:55
Hi,

I am trying to run the spatial empirical orthogonal function on a raster stack using eof function in spacetime 1.2-1 package.

My raster stack is

> stk
class       : RasterStack
dimensions  : 2423, 2470, 5984810, 690  (nrow, ncol, ncell, nlayers)
resolution  : 499.9339, 499.9718  (x, y)
extent      : -85163.44, 1149673, 5656775, 6868207  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : X2001.01.01, X2001.01.09, X2001.01.17, X2001.01.25, X2001.02.02, X2001.02.10, X2001.02.18, X2001.02.26, X2001.03.06, X2001.03.14, X2001.03.22, X2001.03.30, X2001.04.07, X2001.04.15, X2001.04.23, ...
min values  :           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0, ...
max values  :         100,         100,         100,         100,         100,         100,         100,         100,         100,         100,         100,         100,         100,         100,         100, ...

>object.size(stk)
8699416 bytes

The size of raster stack is ~8 mb but when I try to convert to STFDF it become several GBs until I receive error
about memory.

My system is
sysname        release        version        machine
"Windows"     ">= 8 x64"    "build 9200"      "x86-64"

>memory.limit()
[1] 32659

My R version is
"R version 3.5.0 (2018-04-23)"

I would appreciate any advise on possible solutions.

Thank you.

Shiva


        [[alternative HTML version deleted]]

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

spatio temporal anisotropy: help with stAni argument

Mon, 05/21/2018 - 08:16
Hi,
I have some problems in understanding the usage of argument 'stAni' of
vgmST and krigeST functions, for.the setting of spatio temporal anisotropy.
How have I to interpret this argument, that is the integer value that I saw
I have to assign to it, what does it means?
How can this argument influence the result of prediction?
Are units of spatial and temporal ranges and lags related to this argument?
How can I study the spatio temporal anisotropy of my data?
Kind regards.

        [[alternative HTML version deleted]]

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

Re: Intersection of polygons with raster

Sat, 05/19/2018 - 12:40
Thank you Michael, that is exactly what I needed. It is fairly slow for my
case but I only need to do it once :-)

And thank you for providing an example, I'm not too familiar with rasters
and would have had some trouble creating it.

Kent

On Fri, May 18, 2018 at 11:41 PM, Michael Sumner <[hidden email]> wrote:

> Hi Kent, this is pretty straightforward with raster::extract. If speed is
> an issue you can burn the polygon id into the grid and then group values by
> that.
>
> library(raster)
> r <- raster(volcano)
>
> ## very simplistic polygon layer example
> r2 <- raster(r)
> res(r2) <- res(r2) * 20
> p <- sf::st_as_sf(as(r2, "SpatialPolygonsDataFrame"))
>
> ## a dummy threshold value
> onefoot <- 150
> ## grab the first column from extract (it gives multiple columns for
> multi-layer rasters)
> p$onefoot <- extract(r, p, fun = function(x, na.rm = TRUE) any(x >
> onefoot))[,1]
>
> plot(sf::st_geometry(p), col = p$onefoot + 1)
> contour(r, levels = onefoot, col = "white", add = T)
>
>
> ## if speed is an issue, use fasterize for a more abstract workflow
> library(fasterize)
> p$rownum <- 1:nrow(p)
> idgrid <- fasterize(p, r, field = "rownum")
> p$onefoot <- tapply(values(r), values(idgrid), function(x, na.rm = TRUE)
> any(x > onefoot))
> plot(p)
>
> FWIW it's always useful to provide a reprex, since that does take time and
> might miss the mark for your situation. It's also trickier to provide two
> objects that are relevant to each other, so any upfront work you can do in
> an example helps the answerer.
>
> HTH
>
> On Sat, 19 May 2018 at 12:20 Kent Johnson <[hidden email]> wrote:
>
>> Hi,
>>
>> I have a raster object `flood` containing projected flooding levels and a
>> simple features object `parcels` containing MULTIPOLYGONs representing
>> property parcels. I would like to find all the parcel polygons for which
>> there is any flood > 1 foot. What function(s) can do this? Something like
>> raster::intersect(parcels, flood >= 1)
>>
>> I'll put together a reprex if that helps. Hoping someone can point me to
>> the right function to do this, I haven't found anything promising.
>>
>> Thanks,
>> Kent
>>
>>         [[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
> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
> Kingston Tasmania 7050 Australia
> <https://maps.google.com/?q=203+Channel+Highway+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>
>
        [[alternative HTML version deleted]]

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

Re: Intersection of polygons with raster

Fri, 05/18/2018 - 22:41
Hi Kent, this is pretty straightforward with raster::extract. If speed is
an issue you can burn the polygon id into the grid and then group values by
that.

library(raster)
r <- raster(volcano)

## very simplistic polygon layer example
r2 <- raster(r)
res(r2) <- res(r2) * 20
p <- sf::st_as_sf(as(r2, "SpatialPolygonsDataFrame"))

## a dummy threshold value
onefoot <- 150
## grab the first column from extract (it gives multiple columns for
multi-layer rasters)
p$onefoot <- extract(r, p, fun = function(x, na.rm = TRUE) any(x >
onefoot))[,1]

plot(sf::st_geometry(p), col = p$onefoot + 1)
contour(r, levels = onefoot, col = "white", add = T)


## if speed is an issue, use fasterize for a more abstract workflow
library(fasterize)
p$rownum <- 1:nrow(p)
idgrid <- fasterize(p, r, field = "rownum")
p$onefoot <- tapply(values(r), values(idgrid), function(x, na.rm = TRUE)
any(x > onefoot))
plot(p)

FWIW it's always useful to provide a reprex, since that does take time and
might miss the mark for your situation. It's also trickier to provide two
objects that are relevant to each other, so any upfront work you can do in
an example helps the answerer.

HTH

On Sat, 19 May 2018 at 12:20 Kent Johnson <[hidden email]> wrote:

> Hi,
>
> I have a raster object `flood` containing projected flooding levels and a
> simple features object `parcels` containing MULTIPOLYGONs representing
> property parcels. I would like to find all the parcel polygons for which
> there is any flood > 1 foot. What function(s) can do this? Something like
> raster::intersect(parcels, flood >= 1)
>
> I'll put together a reprex if that helps. Hoping someone can point me to
> the right function to do this, I haven't found anything promising.
>
> Thanks,
> Kent
>
>         [[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

Intersection of polygons with raster

Fri, 05/18/2018 - 21:20
Hi,

I have a raster object `flood` containing projected flooding levels and a
simple features object `parcels` containing MULTIPOLYGONs representing
property parcels. I would like to find all the parcel polygons for which
there is any flood > 1 foot. What function(s) can do this? Something like
raster::intersect(parcels, flood >= 1)

I'll put together a reprex if that helps. Hoping someone can point me to
the right function to do this, I haven't found anything promising.

Thanks,
Kent

        [[alternative HTML version deleted]]

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

Re: Raster extract by polygon generating NAs

Thu, 05/17/2018 - 04:45
Thank you Ben, setting na.rm = T did the job.

Best wishes
Joao

On 16 May 2018 at 19:36, Ben Tupper <[hidden email]> wrote:
> Hi,
>
> It's hard to know without any reproducible code, but you will want to pay
> close attention to the value of the na.rm argument to
> raster::extract(Raster,SpatialPolygons)
>
> See ?extract for all the details.
>
> Cheers,
> Ben
>
> On May 16, 2018, at 2:27 PM, João Carreiras <[hidden email]> wrote:
>
> Greetings!
>
> I'm using the extract command (raster package) with a raster layer (x)
> and a spatial polygons dataframe (y). I'm using it to extract the sum
> of raster values by each spatial polygon. However, I'm getting NAs as
> a result for some polygon IDs, which I know isn't true because that
> doesn't happen in ArcGIS.
>
> Does anyone experienced the same issue?
>
> Thanks!
> Joao
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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
>
> Tick Forecasting: https://eco.bigelow.org/
>
>
>
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Raster extract by polygon generating NAs

Wed, 05/16/2018 - 13:36
Hi,

It's hard to know without any reproducible code, but you will want to pay close attention to the value of the na.rm argument to raster::extract(Raster,SpatialPolygons)  

See ?extract for all the details.

Cheers,
Ben

> On May 16, 2018, at 2:27 PM, João Carreiras <[hidden email]> wrote:
>
> Greetings!
>
> I'm using the extract command (raster package) with a raster layer (x)
> and a spatial polygons dataframe (y). I'm using it to extract the sum
> of raster values by each spatial polygon. However, I'm getting NAs as
> a result for some polygon IDs, which I know isn't true because that
> doesn't happen in ArcGIS.
>
> Does anyone experienced the same issue?
>
> Thanks!
> Joao
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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

Tick Forecasting: https://eco.bigelow.org/





        [[alternative HTML version deleted]]

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

Raster extract by polygon generating NAs

Wed, 05/16/2018 - 13:27
Greetings!

I'm using the extract command (raster package) with a raster layer (x)
and a spatial polygons dataframe (y). I'm using it to extract the sum
of raster values by each spatial polygon. However, I'm getting NAs as
a result for some polygon IDs, which I know isn't true because that
doesn't happen in ArcGIS.

Does anyone experienced the same issue?

Thanks!
Joao

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

Re: Will upcoming ggplot2::calc() clash with raster::calc()?

Wed, 05/16/2018 - 13:05
We are tracking at https://github.com/tidyverse/ggplot2/issues/2614 —
and have come up with a name that I prefer to calc(): stat().

Hadley

On Wed, May 16, 2018 at 9:42 AM, Barry Rowlingson
<[hidden email]> wrote:
> It seems the options are:
>
> 1. ggplot2 and raster use calc - scripts will have to use raster::calc or
> ggplot2::calc to be ambiguous.
>
>  This is the painful solution. Scripts will break. Users will have to type
> raster::calc or ggplot2::calc depending on the order they do
> library(raster);library(ggplot2). Both packages are often used together.
>
> 2. `calc` becomes a "generic" and despatches to the right specific function
> when called on a raster or whatever ggplot2::calc needs.
>
>  This either requires "ggplot2" to have "raster" as a dependency or for the
> generic to be moved to a package that each of these will depend on. An S3
> or S4 generic? Will it even work? Nobody wants to manage the dependency.
>
> 3. Use a different name in `ggplot2`, for example, `val`, which has the
> same sense - compare:
>
>    geom_histogram(aes(y = *calc*(count)))
>    geom_histogram(aes(y = *val*(count)))
>
> 4. ggplot2 does something else - maybe the aesthetic could be specified as
> a formula?
>
>   geom_histogram(aes(y = ~count))
>
>  but that would require a chunk of ggplot2 rewriting
>
> Barry
>
>
>
>
>
> On Tue, May 15, 2018 at 6:58 PM, Edzer Pebesma <
> [hidden email]> wrote:
>
>> This was asked on twitter:
>> https://twitter.com/hadleywickham/status/996430251499536385 , by Hadley
>> Wickham:
>>
>> "Can any raster user comment on if the new ggplot2::calc() function is
>> going to cause significant pain because it clashes with raster::calc()?"
>>
>> I'll report answers back to twitter, or directly to Hadley if it doesn't
>> fit in 280 chars.
>> --
>> Edzer Pebesma
>> Institute for Geoinformatics
>> Heisenbergstrasse 2, 48151 Muenster, Germany
>> Phone: +49 251 8333081
>>
>> _______________________________________________
>> 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


--
http://hadley.nz

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

Re: Will upcoming ggplot2::calc() clash with raster::calc()?

Wed, 05/16/2018 - 11:42
It seems the options are:

1. ggplot2 and raster use calc - scripts will have to use raster::calc or
ggplot2::calc to be ambiguous.

 This is the painful solution. Scripts will break. Users will have to type
raster::calc or ggplot2::calc depending on the order they do
library(raster);library(ggplot2). Both packages are often used together.

2. `calc` becomes a "generic" and despatches to the right specific function
when called on a raster or whatever ggplot2::calc needs.

 This either requires "ggplot2" to have "raster" as a dependency or for the
generic to be moved to a package that each of these will depend on. An S3
or S4 generic? Will it even work? Nobody wants to manage the dependency.

3. Use a different name in `ggplot2`, for example, `val`, which has the
same sense - compare:

   geom_histogram(aes(y = *calc*(count)))
   geom_histogram(aes(y = *val*(count)))

4. ggplot2 does something else - maybe the aesthetic could be specified as
a formula?

  geom_histogram(aes(y = ~count))

 but that would require a chunk of ggplot2 rewriting

Barry





On Tue, May 15, 2018 at 6:58 PM, Edzer Pebesma <
[hidden email]> wrote:

> This was asked on twitter:
> https://twitter.com/hadleywickham/status/996430251499536385 , by Hadley
> Wickham:
>
> "Can any raster user comment on if the new ggplot2::calc() function is
> going to cause significant pain because it clashes with raster::calc()?"
>
> I'll report answers back to twitter, or directly to Hadley if it doesn't
> fit in 280 chars.
> --
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
>
> _______________________________________________
> 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

Will upcoming ggplot2::calc() clash with raster::calc()?

Tue, 05/15/2018 - 12:58
This was asked on twitter:
https://twitter.com/hadleywickham/status/996430251499536385 , by Hadley
Wickham:

"Can any raster user comment on if the new ggplot2::calc() function is
going to cause significant pain because it clashes with raster::calc()?"

I'll report answers back to twitter, or directly to Hadley if it doesn't
fit in 280 chars.
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

Resources and Technical issues about Spatio Temporal Kriging

Mon, 05/14/2018 - 08:52
Hi,
I am trying to do spatio temporal kriging. I have produced some prediction
maps, but during the process, I noticed something of strange:
during the execution of the function krigeST I saw a high usage of RAM
memory and Disk, which sometimes produced the crash of RStudio. It is
normal? I am using a spatio-temporal grid of about 2000 points as 'newdata'
parameter in krigeST function.
Another issues is that the result of the function krigeST seems to be a
multiple map (for example I made a temporal grid of length 2, so I got 2 in
1 prediction  map as the result of krigeST). I am not sure if the "number"
of prediction maps depends on the length of temporal grid. In any way, it
could.be possible to obtain a unique prediction map as the result of
krigeST (as in the pure spatial context, for example)?

Kind regards.

        [[alternative HTML version deleted]]

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

Pages