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: 1 hour 40 min ago

Re: Extract coordinates from rasterbrick

Tue, 08/14/2018 - 14:32
Dear Ákos and Vijay,

Thank you very much for the solutions. They work perfectly!

Sincerely,

Milu

On Tue, Aug 14, 2018 at 6:12 PM, Vijay Lulla <[hidden email]> wrote:

> And if you need coordinates as part of your data frame just do
> cbind(coordinates(x), as.data.frame(x, xy = TRUE))
>
> On Tue, Aug 14, 2018 at 11:50 AM Bede-Fazekas Ákos <[hidden email]>
> wrote:
>
> > Dear Milu,
> >
> > I think that you are looking for as.data.frame(x, xy = TRUE).
> >
> > HTH,
> > Ákos Bede-Fazekas
> > Hungarian Academy of Sciences
> >
> >
> > 2018.08.14. 17:03 keltezéssel, Miluji Sb írta:
> > > Dear all,
> > >
> > > I have the following rasterbrick (x);
> > >
> > > class       : RasterBrick
> > > dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
> > > resolution  : 0.25, 0.25  (x, y)
> > > extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
> > > coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > > data source : /projectnb/climpct/TEX_1986_2005.nc
> > > names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
> > > X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
> > > X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
> > > X1986.01.15, ...
> > > Date        : 1986-01-01, 2005-12-31 (min, max)
> > > varname     : tex
> > >
> > >
> > > I can obtain the values by;
> > >
> > > as.data.frame(getValues(x))
> > >
> > > However, I would like to extract the values by all coordinates in the
> > file
> > > in the form of a dataframe. Is it possible to do so? Any help will be
> > > greatly appreciated. Thank you.
> > >
> > > Sincerely,
> > >
> > > Milu
> > >
> > >       [[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
>
        [[alternative HTML version deleted]]

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

Re: Extract coordinates from rasterbrick

Tue, 08/14/2018 - 11:12
And if you need coordinates as part of your data frame just do
cbind(coordinates(x), as.data.frame(x, xy = TRUE))

On Tue, Aug 14, 2018 at 11:50 AM Bede-Fazekas Ákos <[hidden email]>
wrote:

> Dear Milu,
>
> I think that you are looking for as.data.frame(x, xy = TRUE).
>
> HTH,
> Ákos Bede-Fazekas
> Hungarian Academy of Sciences
>
>
> 2018.08.14. 17:03 keltezéssel, Miluji Sb írta:
> > Dear all,
> >
> > I have the following rasterbrick (x);
> >
> > class       : RasterBrick
> > dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
> > resolution  : 0.25, 0.25  (x, y)
> > extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
> > coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > data source : /projectnb/climpct/TEX_1986_2005.nc
> > names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
> > X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
> > X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
> > X1986.01.15, ...
> > Date        : 1986-01-01, 2005-12-31 (min, max)
> > varname     : tex
> >
> >
> > I can obtain the values by;
> >
> > as.data.frame(getValues(x))
> >
> > However, I would like to extract the values by all coordinates in the
> file
> > in the form of a dataframe. Is it possible to do so? Any help will be
> > greatly appreciated. Thank you.
> >
> > Sincerely,
> >
> > Milu
> >
> >       [[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

Re: Extract coordinates from rasterbrick

Tue, 08/14/2018 - 10:49
Dear Milu,

I think that you are looking for as.data.frame(x, xy = TRUE).

HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences


2018.08.14. 17:03 keltezéssel, Miluji Sb írta:
> Dear all,
>
> I have the following rasterbrick (x);
>
> class       : RasterBrick
> dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
> resolution  : 0.25, 0.25  (x, y)
> extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : /projectnb/climpct/TEX_1986_2005.nc
> names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
> X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
> X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
> X1986.01.15, ...
> Date        : 1986-01-01, 2005-12-31 (min, max)
> varname     : tex
>
>
> I can obtain the values by;
>
> as.data.frame(getValues(x))
>
> However, I would like to extract the values by all coordinates in the file
> in the form of a dataframe. Is it possible to do so? Any help will be
> greatly appreciated. Thank you.
>
> Sincerely,
>
> Milu
>
> [[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

Extract coordinates from rasterbrick

Tue, 08/14/2018 - 10:03
Dear all,

I have the following rasterbrick (x);

class       : RasterBrick
dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : /projectnb/climpct/TEX_1986_2005.nc
names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
X1986.01.15, ...
Date        : 1986-01-01, 2005-12-31 (min, max)
varname     : tex


I can obtain the values by;

as.data.frame(getValues(x))

However, I would like to extract the values by all coordinates in the file
in the form of a dataframe. Is it possible to do so? Any help will be
greatly appreciated. Thank you.

Sincerely,

Milu

        [[alternative HTML version deleted]]

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

Re: How do i plot a graph with x axis vertically labelled?

Mon, 08/13/2018 - 09:32
The graphics parameters are all in the help file for par
?par

The one you need is las.

Without trying to wade thru your mangled HTML data (dput() is the best
way to provide reproducible data), check this out:

barplot(runif(5), names.arg = c("25/08/2004", "26/09/2004",
"02/11/2004", "11/11/2004", "22/11/2004"), main = 'Females 2005', ylab
= ' Age/Days', xlab = 'Birth/Hatch date', las=2)

You'll need to mess with label positioning, and with margin size, to
get everything looking as you wish.

But I'm not sure this is the most effective way to show
irregularly-spaced dates. Why not convert your dates to Date or other
time series, and plot them to scale as points?

Sarah
On Mon, Aug 13, 2018 at 6:25 AM Vasana Tutjavi via R-sig-Geo
<[hidden email]> wrote:
>
> Dear R-sig-geo experts,
>
> I want to plot a bar graph with the labels on the x axis in a vertical format. The data is the Ages of Squid (on the Y axis) and the date of Birth on the X axis, I am struggling to change the labels on the x axis to a vertical format and I want all the dates to show at all the bars. My data is as follows:
> Birth
> Sex
> Age
> 25/08/2004
> F
> 2
> 26/09/2004
> F
> 3
> 02/11/2004
> F
> 3
> 11/11/2004
> F
> 1
> 22/11/2004
> F
> 3
> 04/12/2004
> F
> 1
> 08/12/2004
> F
> 4
> 15/12/2004
> F
> 4
> 21/12/2004
> F
> 5
> 24/12/2004
> F
> 1
> 04/01/2005
> F
> 1
> 06/01/2005
> F
> 3
> 13/01/2005
> F
> 2
> 18/01/2005
> F
> 1
> 19/01/2005
> F
> 3
> 09/02/2005
> F
> 2
> 11/02/2005
> F
> 1
> 13/02/2005
> F
> 1
> 17/02/2005
> F
> 1
> 18/02/2005
> F
> 8
> 19/02/2005
> F
> 5
> 24/02/2005
> F
> 1
> 28/02/2005
> F
> 42
> 01/03/2005
> F
> 1
> 07/03/2005
> F
> 1
> 08/03/2005
> F
> 2
> 09/03/2005
> F
> 1
> 14/03/2005
> F
> 7
> 16/03/2005
> F
> 1
> 25/03/2005
> F
> 3
> 01/04/2005
> F
> 1
> 06/04/2005
> F
> 5
> 07/04/2005
> F
> 1
> 08/04/2005
> F
> 1
> 13/04/2005
> F
> 1
> 14/04/2005
> F
> 1
> 27/04/2005
> F
> 1
> 01/05/2005
> F
> 2
> 07/05/2005
> F
> 1
> 08/05/2005
> F
> 2
> 16/05/2005
> F
> 1
> 30/05/2005
> F
> 8
> 06/06/2005
> F
> 1
> Here is the script I tried:
>
> barplot(Females2005$Age,  names.arg = c("25/08/2004", "26/09/2004", "02/11/2004", "11/11/2004", "22/11/2004", "04/12/2004", "08/12/2004",
>                                        "15/12/2004", "21/12/2004", "24/12/2004", "04/01/2005", "06/01/2005", "13/01/2005", "18/01/2005",
>                                        "19/01/2005", "09/02/2005", "11/02/2005", "13/02/2005", "17/02/2005", "18/02/2005", "19/02/2005",
>                                       "24/02/2005", "28/02/2005", "01/03/2005", "07/03/2005", "08/03/2005", "09/03/2005", "14/03/2005",
>                                     "16/03/2005", "25/03/2005", "01/04/2005", "06/04/2005", "07/04/2005", "08/04/2005", "13/04/2005",
>                                       "14/04/2005", "27/04/2005", "01/05/2005", "07/05/2005", "08/05/2005", "16/05/2005", "30/05/2005",
>                                     "06/06/2005"), main = 'Females 2005', ylab = 'Age/Days', xlab = 'Birth/Hatch date')
>
> Many thanks
> Vasana
> Sent from Mail for Windows 10
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


--
Sarah Goslee
http://www.functionaldiversity.org

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

Re: How do i plot a graph with x axis vertically labelled?

Mon, 08/13/2018 - 09:28
Vasana--

May I ask why you want to produce such a graph? What relationship(s) in
your data are you trying to show?  I can't be sure from your original post,
but it seems to me that you are trying to do this the hard way.  R can help
you.

What happens when you run the code you posted? What error message do you
get, if any?

Your birthdates are not just nominal labels; they contain more information
than that--specifically, they are ordered, with certain elapsed duration
between them. Using them merely as names/labels on bars sacrifices
information. Also, if you don't have an observation for every birthdate in
a sequence of dates, then you will represent unequally-spaced dates with
equally-spaced bars--generally not a good thing. Conversely, could any of
your birthdates ever be duplicated?

Vertically-oriented bar labels can be difficult for your intended audience
to read, especially since you seem to want to put 30-40 of them along the
horizontal axis. Another option would be to flip your graph, putting
birthdates on the vertical axis so the date labels can be horizontal.

Here is a minimal working example to demonstrate one possible solution that
you might be able to modify to convey what you are trying to convey, and
that let's R do the work it is good at.

## generate some data
dd <- data.frame( birth = Sys.Date() + sort(sample(1:10, 6, replace =
TRUE)),  sex = factor(sample(c("M", "F"), 6, replace = TRUE)), age =
sample(1:9, 6, replace = TRUE) )
## show the data
dd
str(dd)
## a plot
 with(dd, plot(age ~ birth, type = "h"))

--Chris Ryan

On Mon, Aug 13, 2018 at 6:24 AM, Vasana Tutjavi via R-sig-Geo <
[hidden email]> wrote:

> Dear R-sig-geo experts,
>
> I want to plot a bar graph with the labels on the x axis in a vertical
> format. The data is the Ages of Squid (on the Y axis) and the date of Birth
> on the X axis, I am struggling to change the labels on the x axis to a
> vertical format and I want all the dates to show at all the bars. My data
> is as follows:
> Birth
> Sex
> Age
> 25/08/2004
> F
> 2
> 26/09/2004
> F
> 3
> 02/11/2004
> F
> 3
> 11/11/2004
> F
> 1
> 22/11/2004
> F
> 3
> 04/12/2004
> F
> 1
> 08/12/2004
> F
> 4
> 15/12/2004
> F
> 4
> 21/12/2004
> F
> 5
> 24/12/2004
> F
> 1
> 04/01/2005
> F
> 1
> 06/01/2005
> F
> 3
> 13/01/2005
> F
> 2
> 18/01/2005
> F
> 1
> 19/01/2005
> F
> 3
> 09/02/2005
> F
> 2
> 11/02/2005
> F
> 1
> 13/02/2005
> F
> 1
> 17/02/2005
> F
> 1
> 18/02/2005
> F
> 8
> 19/02/2005
> F
> 5
> 24/02/2005
> F
> 1
> 28/02/2005
> F
> 42
> 01/03/2005
> F
> 1
> 07/03/2005
> F
> 1
> 08/03/2005
> F
> 2
> 09/03/2005
> F
> 1
> 14/03/2005
> F
> 7
> 16/03/2005
> F
> 1
> 25/03/2005
> F
> 3
> 01/04/2005
> F
> 1
> 06/04/2005
> F
> 5
> 07/04/2005
> F
> 1
> 08/04/2005
> F
> 1
> 13/04/2005
> F
> 1
> 14/04/2005
> F
> 1
> 27/04/2005
> F
> 1
> 01/05/2005
> F
> 2
> 07/05/2005
> F
> 1
> 08/05/2005
> F
> 2
> 16/05/2005
> F
> 1
> 30/05/2005
> F
> 8
> 06/06/2005
> F
> 1
> Here is the script I tried:
>
> barplot(Females2005$Age,  names.arg = c("25/08/2004", "26/09/2004",
> "02/11/2004", "11/11/2004", "22/11/2004", "04/12/2004", "08/12/2004",
>                                        "15/12/2004", "21/12/2004",
> "24/12/2004", "04/01/2005", "06/01/2005", "13/01/2005", "18/01/2005",
>                                        "19/01/2005", "09/02/2005",
> "11/02/2005", "13/02/2005", "17/02/2005", "18/02/2005", "19/02/2005",
>                                       "24/02/2005", "28/02/2005",
> "01/03/2005", "07/03/2005", "08/03/2005", "09/03/2005", "14/03/2005",
>                                     "16/03/2005", "25/03/2005",
> "01/04/2005", "06/04/2005", "07/04/2005", "08/04/2005", "13/04/2005",
>                                       "14/04/2005", "27/04/2005",
> "01/05/2005", "07/05/2005", "08/05/2005", "16/05/2005", "30/05/2005",
>                                     "06/06/2005"), main = 'Females 2005',
> ylab = 'Age/Days', xlab = 'Birth/Hatch date')
>
> Many thanks
> Vasana
> Sent from Mail for Windows 10
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

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

Re: [localmoran.exact]

Mon, 08/13/2018 - 09:12
On Mon, 13 Aug 2018, Lisa Menez wrote:

>
> Dear all,
>
> I have been trying to obtain local Moran’s I with the exact approach
> (localmoran.exact) but I have some issues and would need your help.

A complete reproducible example is required (script and data, data may be
downloaded, do not attach). Does the same thing happen with
localmoran.sad()?

>
> I am trying to bring to light how spatial patterns of wealth have
> evolved from 2000 to 2016 in the European Union.
>
> Looking at the Regional Gross Domestic Product (GDP) of European regions, my lm model looks like this :
>
> lm(log(EU_NUTS@data[,year[i]])~1+factor(EU_NUTS@data$NUTS_L1))
>
> where log(EU_NUTS@data[,year[i]]) are log of regional GDPs
> Never, ever use @ to access data. EU_NUTS@data$NUTS_L1 should be
EU_NUTS$NUTS_L1 and should be (made a) factor first. It isn't clear what
log(EU_NUTS@data[,year[i]]) will end up being - is "year" a character
vector? If so, EU_NUTS[[year[i]]] is the correct syntax.

> and factor(EU_NUTS@data$NUTS_L1) is the vector of country dummies.
>
> Sadly, I obtain an error message that says :
>
> Error in integrate(integrand, lower = 0, upper = upper) :
>  a limit is missing
> Moreover : Warning messages:
> 1: In sqrt(2 * t2 - t1^2) : production of NaN
> 2: In sqrt(2 * t2 - t1^2) : production of NaN
>
> According to the GitHub source code , No, as you can see from https://cran.r-project.org/package=spdep, the
development site of spdep is https://github.com/r-spatial/spdep/. Never
use non-authorised copies. To see the source, simply type the function
name. Use debug(localmoran.exact) to step through the function while it
runs, to see the values of the objects you are in doubt about. Do not do
this in RStudio, its debugging interface tries to be far too clever.

The following is illegible in plain text; post in plain text only.

Roger

>
>        Vi <- listw2star <https://rdrr.io/cran/spdep/man/localmoran.sad.html>(B, select[i], style=style, n, D <https://rdrr.io/r/stats/deriv.html>, a,
>    zero.policy=zero.policy)
>        Viu <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, u, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
> Ii <- c <https://rdrr.io/r/base/c.html>(crossprod <https://rdrr.io/r/base/crossprod.html>(u, Viu) / utu)
>        ViX <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, X, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
>        MViM <- t <https://rdrr.io/r/base/t.html>(X) %*% ViX %*% XtXinv
>        t1 <- -sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM))
>        sumsq.Vi <- function <https://rdrr.io/r/base/function.html>(x) {
>            if <https://rdrr.io/r/base/Control.html> (is.null <https://rdrr.io/r/base/NULL.html>(x)) NA <https://rdrr.io/r/base/NA.html>
>    else <https://rdrr.io/r/base/Control.html> sum <https://rdrr.io/r/base/sum.html>(x^2)
> }
> trVi2 <- sum <https://rdrr.io/r/base/sum.html>(sapply <https://rdrr.io/r/base/lapply.html>(Vi$weights <https://rdrr.io/r/stats/weights.html>, sumsq.Vi), na.rm=TRUE <https://rdrr.io/r/base/logical.html>)
> t2a <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(t <https://rdrr.io/r/base/t.html>(ViX) %*% ViX %*% XtXinv))
> t2b <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM %*% MViM))
> t2 <- trVi2 - 2*t2a + t2b
> e1 <- 0.5 * (t1 + sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
> en <- 0.5 * (t1 - sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
>
> My guess is that my setting causes (2*t2<t1^2) : Is it ? Why would it ?
>
> My database is available here : http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt>
> and my code : http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt>
>
> I thank you very much for your attention and your help.
>
> Best Regards
> Lisa
>
> PhD Student
> GREDEG, Université Côté d’Azur
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

[localmoran.exact]

Mon, 08/13/2018 - 08:56

Dear all,

I have been trying to obtain local Moran’s I with the exact approach (localmoran.exact)
 but I have some issues and would need your help.

I am trying to bring to light how spatial patterns of wealth have evolved from 2000 to 2016 in the European Union.

Looking at the Regional Gross Domestic Product (GDP) of European regions, my lm model looks like this :

lm(log(EU_NUTS@data[,year[i]])~1+factor(EU_NUTS@data$NUTS_L1))

where log(EU_NUTS@data[,year[i]]) are log of regional GDPs

and factor(EU_NUTS@data$NUTS_L1) is the vector of country dummies.

Sadly, I obtain an error message that says :

Error in integrate(integrand, lower = 0, upper = upper) :
  a limit is missing
Moreover : Warning messages:
1: In sqrt(2 * t2 - t1^2) : production of NaN
2: In sqrt(2 * t2 - t1^2) : production of NaN

According to the GitHub source code ,

        Vi <- listw2star <https://rdrr.io/cran/spdep/man/localmoran.sad.html>(B, select[i], style=style, n, D <https://rdrr.io/r/stats/deriv.html>, a,
            zero.policy=zero.policy)
        Viu <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, u, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
        Ii <- c <https://rdrr.io/r/base/c.html>(crossprod <https://rdrr.io/r/base/crossprod.html>(u, Viu) / utu)
        ViX <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, X, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
        MViM <- t <https://rdrr.io/r/base/t.html>(X) %*% ViX %*% XtXinv
        t1 <- -sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM))
        sumsq.Vi <- function <https://rdrr.io/r/base/function.html>(x) {
            if <https://rdrr.io/r/base/Control.html> (is.null <https://rdrr.io/r/base/NULL.html>(x)) NA <https://rdrr.io/r/base/NA.html>
            else <https://rdrr.io/r/base/Control.html> sum <https://rdrr.io/r/base/sum.html>(x^2)
        }
        trVi2 <- sum <https://rdrr.io/r/base/sum.html>(sapply <https://rdrr.io/r/base/lapply.html>(Vi$weights <https://rdrr.io/r/stats/weights.html>, sumsq.Vi), na.rm=TRUE <https://rdrr.io/r/base/logical.html>)
        t2a <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(t <https://rdrr.io/r/base/t.html>(ViX) %*% ViX %*% XtXinv))
        t2b <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM %*% MViM))
        t2 <- trVi2 - 2*t2a + t2b
        e1 <- 0.5 * (t1 + sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
        en <- 0.5 * (t1 - sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))

My guess is that my setting causes (2*t2<t1^2) : Is it ? Why would it ?

My database is available here : http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt>
and my code : http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt>

I thank you very much for your attention and your help.

Best Regards
Lisa

PhD Student
GREDEG, Université Côté d’Azur






        [[alternative HTML version deleted]]

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

Question: How to parametrize rstoolbox for radCor operations of ASTERLB1 data

Mon, 08/13/2018 - 05:45
 
Hello good morning learned community of specialists. I would be glad to receive some guidance or code snippets to help me carry out radCor operations on ASTERLB1 data.
Thank you so much.    On Monday, August 13, 2018, 6:06:32 PM GMT+8, [hidden email] <[hidden email]> wrote:  
 
 Send R-sig-Geo mailing list submissions to
    [hidden email]

To subscribe or unsubscribe via the World Wide Web, visit
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
    [hidden email]

You can reach the person managing the list at
    [hidden email]

When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-Geo digest..."


Today's Topics:

  1. Autocorrelation and stationarity in spatio temporal data
      (Bruno Sesti)

----------------------------------------------------------------------

Message: 1
Date: Sun, 12 Aug 2018 17:45:31 +0200
From: Bruno Sesti <[hidden email]>
To: [hidden email]
Subject: [R-sig-Geo] Autocorrelation and stationarity in spatio
    temporal data
Message-ID:
    <[hidden email]>
Content-Type: text/plain; charset="utf-8"

Hi, I am working with spatio temoral data and I am experimenting spatio
temporal data analysis and interpolation. I would like to ask if there is
some technique/package/function in R to check if irregilar spatio temporal
data, in particolar air quality spatio temporal data from mobile sensors,
are stationary (second order stationary or intrinsic stationary) and if the
shows some form of (auto)correlation.

Kind regards.

    [[alternative HTML version deleted]]




------------------------------

Subject: Digest Footer

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


------------------------------

End of R-sig-Geo Digest, Vol 180, Issue 10
******************************************
 
        [[alternative HTML version deleted]]

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

How do i plot a graph with x axis vertically labelled?

Mon, 08/13/2018 - 05:24
Dear R-sig-geo experts, 

I want to plot a bar graph with the labels on the x axis in a vertical format. The data is the Ages of Squid (on the Y axis) and the date of Birth on the X axis, I am struggling to change the labels on the x axis to a vertical format and I want all the dates to show at all the bars. My data is as follows:
Birth
Sex
Age
25/08/2004
F
2
26/09/2004
F
3
02/11/2004
F
3
11/11/2004
F
1
22/11/2004
F
3
04/12/2004
F
1
08/12/2004
F
4
15/12/2004
F
4
21/12/2004
F
5
24/12/2004
F
1
04/01/2005
F
1
06/01/2005
F
3
13/01/2005
F
2
18/01/2005
F
1
19/01/2005
F
3
09/02/2005
F
2
11/02/2005
F
1
13/02/2005
F
1
17/02/2005
F
1
18/02/2005
F
8
19/02/2005
F
5
24/02/2005
F
1
28/02/2005
F
42
01/03/2005
F
1
07/03/2005
F
1
08/03/2005
F
2
09/03/2005
F
1
14/03/2005
F
7
16/03/2005
F
1
25/03/2005
F
3
01/04/2005
F
1
06/04/2005
F
5
07/04/2005
F
1
08/04/2005
F
1
13/04/2005
F
1
14/04/2005
F
1
27/04/2005
F
1
01/05/2005
F
2
07/05/2005
F
1
08/05/2005
F
2
16/05/2005
F
1
30/05/2005
F
8
06/06/2005
F
1
Here is the script I tried:

barplot(Females2005$Age,  names.arg = c("25/08/2004", "26/09/2004", "02/11/2004", "11/11/2004", "22/11/2004", "04/12/2004", "08/12/2004",
                                       "15/12/2004", "21/12/2004", "24/12/2004", "04/01/2005", "06/01/2005", "13/01/2005", "18/01/2005",
                                       "19/01/2005", "09/02/2005", "11/02/2005", "13/02/2005", "17/02/2005", "18/02/2005", "19/02/2005",
                                      "24/02/2005", "28/02/2005", "01/03/2005", "07/03/2005", "08/03/2005", "09/03/2005", "14/03/2005",
                                    "16/03/2005", "25/03/2005", "01/04/2005", "06/04/2005", "07/04/2005", "08/04/2005", "13/04/2005",
                                      "14/04/2005", "27/04/2005", "01/05/2005", "07/05/2005", "08/05/2005", "16/05/2005", "30/05/2005",
                                    "06/06/2005"), main = 'Females 2005', ylab = 'Age/Days', xlab = 'Birth/Hatch date')

Many thanks
Vasana
Sent from Mail for Windows 10


        [[alternative HTML version deleted]]

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

Autocorrelation and stationarity in spatio temporal data

Sun, 08/12/2018 - 10:45
Hi, I am working with spatio temoral data and I am experimenting spatio
temporal data analysis and interpolation. I would like to ask if there is
some technique/package/function in R to check if irregilar spatio temporal
data, in particolar air quality spatio temporal data from mobile sensors,
are stationary (second order stationary or intrinsic stationary) and if the
shows some form of (auto)correlation.

Kind regards.

        [[alternative HTML version deleted]]

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

Re: rgdal compile confusion

Sat, 08/11/2018 - 09:38
On Sat, 11 Aug 2018, Roger Bivand wrote:

> The gcc you have is ancient; I would guess that your GDAL was not built with
> it.

Roger,

   On Slackware everything is built from source other than a few applications
that are repacked by the SlackBuild.org maintainers from .deb or .rpm
packages.

> What happens if you install GDAL from source? Very likely the decision in
> GDAL to go with CXX11 was premature, but tardy platforms need forcing too
> (typical problem on CentOS and RHEL too).

   Slackware stays back from the bleeding edge and upgrades slowly so the
distribution remains stable.

> Could you post on R-sig-geo - there may be someone there in your
> situation.

   Moved to the mail list with this response.

Thanks,

Rich

---------- Original message -----------
>>  Please share your thoughts on what might be going on here.
>>
>> * installing *source* package ‘rgdal’ ...
>> ** package ‘rgdal’ successfully unpacked and MD5 sums checked
>> configure:  R_HOME: /usr/lib/R
>> configure:  CC: gcc
>> configure:  CXX: g++
>> configure:  C++11 support available
>> configure:  rgdal: 1.3-4
>> checking for /usr/bin/svnversion... yes
>> configure: svn revision: 766
>> checking for gdal-config... /usr/bin/gdal-config
>> checking gdal-config usability... yes
>> configure: GDAL: 2.3.0
>> checking C++11 support for GDAL >= 2.3.0... yes
>> checking GDAL version >= 1.11.4... yes
>> checking gdal: linking with --libs only... no
>> checking gdal: linking with --libs and --dep-libs... no
>> In file included from /usr/include/gdal.h:45:0,
>>                 from gdal_test.cc:1:
>> /usr/include/cpl_port.h:187:6: error: #error Must have C++11 or newer.
>> #    error Must have C++11 or newer.
>>      ^
>> In file included from /usr/include/gdal.h:49:0,
>>                 from gdal_test.cc:1:
>> /usr/include/cpl_minixml.h:202:47: error: expected template-name before '<'
>> token
>> class CPLXMLTreeCloser: public std::unique_ptr<CPLXMLNode,
>> CPLXMLTreeCloserDeleter>
>>                                               ^
>> /usr/include/cpl_minixml.h:202:47: error: expected '{' before '<' token
>> /usr/include/cpl_minixml.h:202:47: error: expected unqualified-id before
>> '<' token
>> In file included from /usr/include/ogr_api.h:45:0,
>>                  from /usr/include/gdal.h:50,
>>                  from gdal_test.cc:1:
>> /usr/include/ogr_core.h:79:28: error: expected '}' before end of line
>> /usr/include/ogr_core.h:79:28: error: expected declaration before end of
>> line
>> In file included from /usr/include/gdal.h:45:0,
>>                 from gdal_test.cc:1:
>> /usr/include/cpl_port.h:187:6: error: #error Must have C++11 or newer.
>> #    error Must have C++11 or newer.
>>      ^
>> In file included from /usr/include/gdal.h:49:0,
>>                 from gdal_test.cc:1:
>> /usr/include/cpl_minixml.h:202:47: error: expected template-name before '<'
>> token
>> class CPLXMLTreeCloser: public std::unique_ptr<CPLXMLNode,
>> CPLXMLTreeCloserDeleter>
>>                                               ^
>> /usr/include/cpl_minixml.h:202:47: error: expected '{' before '<' token
>> /usr/include/cpl_minixml.h:202:47: error: expected unqualified-id before
>> '<' token
>> In file included from /usr/include/ogr_api.h:45:0,
>>                  from /usr/include/gdal.h:50,
>>                  from gdal_test.cc:1:
>> /usr/include/ogr_core.h:79:28: error: expected '}' before end of line
>> /usr/include/ogr_core.h:79:28: error: expected declaration before end of
>> line
>> configure:  Install failure: compilation and/or linkage problems.
>> configure:  error: GDALAllRegister not found in libgdal.
>> ERROR: configuration failed for package ‘rgdal’
>> * removing ‘/usr/lib/R/library/rgdal’
>> * restoring previous ‘/usr/lib/R/library/rgdal’
>>
>> The downloaded source packages are in
>> ‘/tmp/RtmpVGIz3Q/downloaded_packages’
>> Updating HTML index of packages in '.Library'
>> Making 'packages.html' ... done
>> Warning message:
>> In install.packages("rgdal") :
>>   installation of package ‘rgdal’ had non-zero exit status
>>
>>  Installed here are:
>> gcc-5.5.0-i586-1_slack14.2
>> gcc-g++-5.5.0-i586-1_slack14.2
>> gcc-gfortran-5.5.0-i586-1_slack14.2
>> gcc-gnat-5.5.0-i586-1_slack14.2
>> gcc-go-5.5.0-i586-1_slack14.2
>> gcc-java-5.5.0-i586-1_slack14.2
>> gcc-objc-5.5.0-i586-1_slack14.2
>> gccmakedep-1.0.3-noarch-1
>>
>> and
>>
>> gdal-2.3.0-i586-1_SBo
>>
>> Regards,
>>
>> Rich
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

www.asdar-book.org site issues

Sat, 08/11/2018 - 06:36
The ASDAR book site is not being shown for viewing, but the content is
still there. Anyone needing content (several inquiries off-list) may try:

ASDAR_BOOK <- "http://www.asdar-book.org/bundles2ed"
chapters <- c("hello", "cm", "vis", "die", "cm2", "std", "sppa", "geos",
  "lat", "dismap")
for (i in chapters) {fn <- paste(i, "bundle.zip", sep="_");
download.file(paste(ASDAR_BOOK, fn, sep = "/"), fn)}

to download the chapter bundles. Summer maintenance takes time, and we
hope to fix the underlying problem in a week or so.

Roger

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

help regarding methodology

Fri, 08/10/2018 - 13:09
I want to work on conditional probability using spatial analysis. To be
more precise,  as in one event depends on the other but not vice versa and
I want to prove this relation spatially. what methodology I can use.

        [[alternative HTML version deleted]]

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

Re: raster::overlay with a function that calls a list of vectors

Thu, 08/09/2018 - 08:22
Hi Rafael.

Your suggestion works perfectly for the toy data.
For my real data, it's not that simple. I'm processing data for the entire
Europe, and thus I split data in tiles. Apparently, Vectorize does not work
when a tile includes just some of the strata (raster r2 in the toy data has
2 strata, but my real data has 7). When tiles include all my strata,
Vectorize seems to work, although slowly (it is processing while writing
this message). Anyway, your message was important to spot this behaviour, I
might need to include some workaround to process the tiles with less
strata, while trying to speed up the processing.

Thank you very much
Hugo

Rafael Wüest <[hidden email]> escreveu no dia quinta, 9/08/2018
à(s) 14:33:

> Hi Hugo
>
> Vectorize() sometimes solves this issue.
>
> foo <- function(x,y,...) {return(rbinom(n=1, size=1, prob=probs[[y]][x]))}
> fooV <- Vectorize(foo)
>
> Using your toy example, check carefully if
>
> o <- overlay(r1,r2,fun = fooV)
>
> gives you what you are looking for.
>
> HTH
> Rafael
>
>
> > On 9 Aug 2018, at 11:51, Hugo Costa <[hidden email]> wrote:
> >
> > Dear list
> >
> > I wonder if it’s possible to use overlay function of raster package with
> a
> > function that calls a list of vectors to perform some calculations based
> on
> > two rasters. So far I just saw examples of functions performing some
> raster
> > algebra without calling external data.
> >
> > Hereafter I provide some toy code to illustrate what I’m trying to do,
> but
> > can also provide some context on my real problem. Specifically, I need to
> > classify each pixel as either zero (absence) or one (presence) of
> housings.
> > The likelihood of housing presence is related to the percentage of
> built-up
> > area covering the pixel (raster ‘r1’ below), and the land cover type
> > (raster ‘r2’ below). This likelihood is known based on reference data,
> > which is stored in a list like ‘probs’ below.
> >
> > library(raster)
> > # continuous and categorical maps
> > r1<-r2<-raster()
> > r1[]<-round(runif(ncell(r1))*100)
> > r2[]<-1
> > r2[1:30000]<-2
> > # probability of housing presence in each stratum
> > prob1<-1:100/100
> > prob2<-log(1:100)/max(log(1:100))
> > # list of probabilities to be used in overlay
> > probs<-list(prob1,prob2)
> > # overlay - not working
> > o<-overlay(r1,r2,fun=function(x,y,...){return(rbinom(n=1, size=1,
> > prob=probs[[y]][x]))})
> >
> > the error is
> > cannot use this formula, probably because it is not vectorized
> >
> > Alternatively to the toy code above, I thought to process each
> categorical
> > class separately and use function calc rather than function overlay (see
> > below). However, this is extremely slow (if not impossible) for large
> > rasters, so I though overlay would be better.
> >
> > # alternative: loop across categorical classes (extremely slow for
> > large rasters)
> > r<-list()for(i in 1:2){
> >  stratum<-r2
> >  stratum[Which(stratum !=i)]<-NA
> >  r[[i]]<-calc(r1, fun=function(x,...){return(rbinom(n=1, size=1,
> > prob=probs[[i]][x]))})
> >  r[[i]]<-mask(r[[i]],stratum)}
> >
> > r<-stack(r)
> > r<-sum(r,na.rm=T)
> >
> > par(mfrow=c(1,3))
> > plot(r1)
> > plot(r2)
> > plot(r)
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> Rafael Wüest Karpati
> http://www.rowueest.net
>
>
        [[alternative HTML version deleted]]

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

Re: raster::overlay with a function that calls a list of vectors

Thu, 08/09/2018 - 07:33
Hi Hugo

Vectorize() sometimes solves this issue.

foo <- function(x,y,...) {return(rbinom(n=1, size=1, prob=probs[[y]][x]))}
fooV <- Vectorize(foo)

Using your toy example, check carefully if

o <- overlay(r1,r2,fun = fooV)

gives you what you are looking for.

HTH
Rafael


> On 9 Aug 2018, at 11:51, Hugo Costa <[hidden email]> wrote:
>
> Dear list
>
> I wonder if it’s possible to use overlay function of raster package with a
> function that calls a list of vectors to perform some calculations based on
> two rasters. So far I just saw examples of functions performing some raster
> algebra without calling external data.
>
> Hereafter I provide some toy code to illustrate what I’m trying to do, but
> can also provide some context on my real problem. Specifically, I need to
> classify each pixel as either zero (absence) or one (presence) of housings.
> The likelihood of housing presence is related to the percentage of built-up
> area covering the pixel (raster ‘r1’ below), and the land cover type
> (raster ‘r2’ below). This likelihood is known based on reference data,
> which is stored in a list like ‘probs’ below.
>
> library(raster)
> # continuous and categorical maps
> r1<-r2<-raster()
> r1[]<-round(runif(ncell(r1))*100)
> r2[]<-1
> r2[1:30000]<-2
> # probability of housing presence in each stratum
> prob1<-1:100/100
> prob2<-log(1:100)/max(log(1:100))
> # list of probabilities to be used in overlay
> probs<-list(prob1,prob2)
> # overlay - not working
> o<-overlay(r1,r2,fun=function(x,y,...){return(rbinom(n=1, size=1,
> prob=probs[[y]][x]))})
>
> the error is
> cannot use this formula, probably because it is not vectorized
>
> Alternatively to the toy code above, I thought to process each categorical
> class separately and use function calc rather than function overlay (see
> below). However, this is extremely slow (if not impossible) for large
> rasters, so I though overlay would be better.
>
> # alternative: loop across categorical classes (extremely slow for
> large rasters)
> r<-list()for(i in 1:2){
>  stratum<-r2
>  stratum[Which(stratum !=i)]<-NA
>  r[[i]]<-calc(r1, fun=function(x,...){return(rbinom(n=1, size=1,
> prob=probs[[i]][x]))})
>  r[[i]]<-mask(r[[i]],stratum)}
>
> r<-stack(r)
> r<-sum(r,na.rm=T)
>
> par(mfrow=c(1,3))
> plot(r1)
> plot(r2)
> plot(r)
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


--
Rafael Wüest Karpati
http://www.rowueest.net

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

model predictions for spatially corrected logistic model

Thu, 08/09/2018 - 04:56
Hi

Does any one have an idea how you make a model prediction for a model that is corrected for spatial autocorrelation? I used the cote-corrected covariate mean in the formation of new data. Is it the normal way you predict a model. I tried doing this but a lot my data gets removed and with the little data it is not able to plot.

Some advise or suggestions would be great.

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

raster::overlay with a function that calls a list of vectors

Thu, 08/09/2018 - 04:51
Dear list

I wonder if it’s possible to use overlay function of raster package with a
function that calls a list of vectors to perform some calculations based on
two rasters. So far I just saw examples of functions performing some raster
algebra without calling external data.

Hereafter I provide some toy code to illustrate what I’m trying to do, but
can also provide some context on my real problem. Specifically, I need to
classify each pixel as either zero (absence) or one (presence) of housings.
The likelihood of housing presence is related to the percentage of built-up
area covering the pixel (raster ‘r1’ below), and the land cover type
(raster ‘r2’ below). This likelihood is known based on reference data,
which is stored in a list like ‘probs’ below.

library(raster)
# continuous and categorical maps
r1<-r2<-raster()
r1[]<-round(runif(ncell(r1))*100)
r2[]<-1
r2[1:30000]<-2
# probability of housing presence in each stratum
prob1<-1:100/100
prob2<-log(1:100)/max(log(1:100))
# list of probabilities to be used in overlay
probs<-list(prob1,prob2)
# overlay - not working
o<-overlay(r1,r2,fun=function(x,y,...){return(rbinom(n=1, size=1,
prob=probs[[y]][x]))})

the error is
cannot use this formula, probably because it is not vectorized

Alternatively to the toy code above, I thought to process each categorical
class separately and use function calc rather than function overlay (see
below). However, this is extremely slow (if not impossible) for large
rasters, so I though overlay would be better.

# alternative: loop across categorical classes (extremely slow for
large rasters)
r<-list()for(i in 1:2){
  stratum<-r2
  stratum[Which(stratum !=i)]<-NA
  r[[i]]<-calc(r1, fun=function(x,...){return(rbinom(n=1, size=1,
prob=probs[[i]][x]))})
  r[[i]]<-mask(r[[i]],stratum)}

r<-stack(r)
r<-sum(r,na.rm=T)

par(mfrow=c(1,3))
plot(r1)
plot(r2)
plot(r)

        [[alternative HTML version deleted]]

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

Re: bottom area of a lake calculation

Wed, 08/08/2018 - 09:49
Hello Sophie,
Yes, what you show was my understanding. I'm sure what you want to do can be done in R, but I'm not sure what packages to use. For doing such a thing, I use GRASS GIS. Regardless, I think you should re-project your data from latitude-longitude to Albers Equal Area or Lambert Conic Conformal, etc., that is area preserving. Eventually, you need to do this to get an area anyhow.
An additional step would be to reproject, then spatially interpolate the point locations you have in your spreadsheet into a grid to Albers Equal Area or Lambert Conic Conformal, etc.
I could take a stab at this and show you...
Regards,Tom

On Wed, Aug 8, 2018 at 10:39 AM, <[hidden email]> wrote:
Dear Tom,

thanks for your reply.
Just to be sure I have explain correctly.
With marmap I can get the green area but I also want to calculate the red one.
I have attached the kind of data I have.


I do not have any GIS software.
I wanted to convert my spherical coordinate to cartesian ones with

x = radius * cos(latitude) * cos(longitude)
y = radius * cos(latitude) * sin(longitude)
z = radius * sin(latitude)
radius= earth_radius-depth : not fully sure...
Do you think it is correct?

The slope is not the same everywhere around the lake.
Is it possible to calculate area with your first explanations?
what R package should I use?

thanks,
Sophie Leblanc







http://sophieleblanc6.wix.com/photographie De : "Thomas Adams"
A : [hidden email],"[hidden email]"
Envoyé: mercredi 8 août 2018 16:07
Objet : Re: [R-sig-Geo] bottom area of a lake calculation
  Sophie,   My apologies, I was thinking in terms of using GRASS GIS... but the idea is the same   Tom   On Wed, Aug 8, 2018 at 10:06 AM, Thomas Adams <[hidden email]> wrote: Sophie,   I think a reasonable approximation might be obtained:   (1) if you were to calculate the slope at each grid point using r.slope.aspect (2) use r.mapcalc to calculate A/cos(slope), where A is the pixel projected area (e.g., 100mx100m = 10000 m^2) -- lat-long coordinates will not work.   Consequently, if the slope were approximately 0, cos(slope) approaches 1 and A/cos(slope) = A; as slope approaches 90 deg, cos(slope) becomes very small, so A becomes very large. This works where the slope direction is orthogonal to a pixel edge, so a correction to account for other slope directions may be needed. But, with smaller grids, the error may be negligible for you...   Tom     On Wed, Aug 8, 2018 at 4:55 AM, <[hidden email]> wrote: <span style="font-family:arial,helvetica,sans-serif; font-size:12px">‌</span>Dear all,<br>
<br>
I would like to calculate the bottom area of a lake. I have already used marmap to calculate projected area, i.e. surface area.<br>
But i have not found anything to calculate a non projected area. I have the coordinates of points (longitude, latitude and depth) I can convert them into xyz. But then I don't know how to calculate the surface I have found polygonal area calculation but I have a 3D object, a polyheron and I think I can not use polygonal area.<br>
Do you have any idea? do you know an R package that could help?<br>
<br>
thanks,<br>
Sophie Leblanc<br>
<br>
<!-- +Signature --><font size="2"><em><a href="http://sophieleblanc6.wix.com/photographie">http://sophieleblanc6.wix.com/photographie</a></em></font> <!-- -Signature -->
        [[alternative HTML version deleted]]

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


   

     


--
Thomas E Adams, III1724 Sage LaneBlacksburg, VA 24060[hidden email] (personal)[hidden email] (work)

1 (513) 739-9512 (cell)

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

Re: bottom area of a lake calculation

Wed, 08/08/2018 - 09:07
Sophie,

My apologies, I was thinking in terms of using GRASS GIS... but the idea is
the same

Tom

On Wed, Aug 8, 2018 at 10:06 AM, Thomas Adams <[hidden email]> wrote:

> Sophie,
>
> I think a reasonable approximation might be obtained:
>
> (1) if you were to calculate the slope at each grid point using
> r.slope.aspect
> (2) use r.mapcalc to calculate A/cos(slope), where A is the pixel
> projected area (e.g., 100mx100m = 10000 m^2) -- lat-long coordinates will
> not work.
>
> Consequently, if the slope were approximately 0, cos(slope) approaches 1
> and A/cos(slope) = A; as slope approaches 90 deg, cos(slope) becomes very
> small, so A becomes very large. This works where the slope direction is
> orthogonal to a pixel edge, so a correction to account for other slope
> directions may be needed. But, with smaller grids, the error may be
> negligible for you...
>
> Tom
>
>
> On Wed, Aug 8, 2018 at 4:55 AM, <[hidden email]> wrote:
>
>> <span style="font-family:arial,helvetica,sans-serif;
>> font-size:12px">‌</span>Dear all,<br>
>> <br>
>> I would like to calculate the bottom area of a lake. I have already used
>> marmap to calculate projected area, i.e. surface area.<br>
>> But i have not found anything to calculate a non projected area. I have
>> the coordinates of points (longitude, latitude and depth) I can convert
>> them into xyz. But then I don't know how to calculate the surface I have
>> found polygonal area calculation but I have a 3D object, a polyheron and I
>> think I can not use polygonal area.<br>
>> Do you have any idea? do you know an R package that could help?<br>
>> <br>
>> thanks,<br>
>> Sophie Leblanc<br>
>> <br>
>> <!-- +Signature --><font size="2"><em><a href="http://sophieleblanc6.wi
>> x.com/photographie">http://sophieleblanc6.wix.com/photographie</a></em></font>
>> <!-- -Signature -->
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
>
>
>
        [[alternative HTML version deleted]]

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

Pages