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 43 min ago

Re: Common Factors test in R

Fri, 04/05/2019 - 14:35
On Fri, 5 Apr 2019, Hulényi Martin wrote:

> Hello,
> I would like to perform a common factor test, conducted using likelihood
> ratio test on spatial error model and spatial durbin model (both in
> panel format).
> I have not found a command in R, that would help me to conduct the test.
> Hence, I am trying to perform the test manuallly using the splm package
> and data available in the splm and plm packages.

> Here is my code:
> library(splm)
> library(plm)
>
> data(Produc, package="plm")
> data(usaww)
>
> Produc <- pdata.frame(Produc, index = c("state", "year"))
>
> usaww<- mat2listw(usaww, style="W")
>
> Produc$slagUnemp <- slag(Produc$unemp, listw = usaww)
>
> durbin <- spml(gsp~unemp + slagUnemp,
>                data=Produc, listw=usaww,  effect = "twoways",
>                 model="within", lag=TRUE, spatial.error = "none",
>                  quiet = FALSE)
> spem  <- spml(gsp ~ unemp,
>                data=Produc, listw=usaww,  effect = "twoways",
>                model="within", lag=FALSE, spatial.error = "b",
>                quiet = FALSE)
> Is it correct to take the last value of the function from the console
> output to compute the likelihood ratio? Without checking whether the likelihoods are compatible (here probably
are), you will not see whether the fitting functions concentrate them,
possibly differently. Here both are lag=TRUE but one has
spatial.error="none", the other "b", so without reading the code, you
can't tell. It would be a good idea if these models had logLik() methods,
because then lmtest::lrtest() should work. It then handles the degrees of
freedom accounting between the models.

Hope this helps,

Roger

> Meaning, in this example, to calculate 2*(function(non-nested) -
> function(nested)) = 2*(10261.79 - 10255.74) = 12.1?
> If it is correct, how can I compute the significance values?
> If it is not correct, is there a better way to compute this?
>
> Autorom tejto spr???vy elektronickej po???ty je Martin Hul???nyi. T???to spr???va je ur???en??? v???lu???ne jej adres???tovi. Inform???cie a ???daje, ktor??? s??? v nej uveden???, alebo ktor??? s??? obsiahnut??? v jej prilo???en???ch s???boroch, m???u by??? inform???ciami alebo ???dajmi chr???nen???mi pod???a platn???ch pr???vnych predpisov v Slovenskej republike. V pr???pade, ak nie ste ur???en??? ako prij???mate??? tejto spr???vy alebo jeho opr???vnen??? z???stupca, upozor???ujeme V???s, ???e inform???cie a ???daje v nej uveden??? nie ste opr???vnen??? sprac???va???, ani ich spr???stupni??? alebo poskytn?????? tretej osobe alebo ich zverejni???. Ak ste nedopatren???m prijali alebo zachytili tuto spr???vu elektronickej po???ty, dovo???ujeme si V???s po???iada???, aby ste ju zaslali sp??? na elektronick??? adresu jej odosielate???a, a aby ste ju n???sledne zmazali zo svojho po??????ta???a a z Va???ej schr???nky elektronickej po???ty. Odosielate??? tejto spr???vy nenesie zodpovednos??? za ???kody sp???soben??? nespr???vnym pou???it???m tejto spr???vy elektronickej po???ty a jej pr???loh.
>
> The author of this e-mail message is Martin Hul???nyi. This message is intended solely for the recipient. The information and data contained therein or contained in its enclosed files may be information or data protected under applicable laws in the Slovak Republic. If you are not designated as the recipient of this message or its authorized representative, we would like to inform you that the information and data contained therein are not authorized to process or make them available to third parties or to disclose them. If you have received or downloaded this e-mail message accidentally, we ask you to send it back to the e-mail address of the sender and then delete it from your computer and from your e-mail. The sender of this message is not responsible for damages caused by incorrect use of this e-mail message and its attachments.
> [eco.jpg]       Pred vytla???en???m tohto mailu zv???te pros???m vplyv na ???ivotn??? prostredie. ???akujeme.
> Please consider the environment before printing this e-mail. Thanks
>
> [[alternative HTML version deleted]]
>
> --
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]
https://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

Re: Creating a color shade between upper and lower bounds of a line plot in R

Fri, 04/05/2019 - 12:40
Ah, I see.

But let's say I had 10 lines on the same plot and wanted to shade green between them. Would this still work, or it only works for 2 lines?
Also, for the "rev" function in there, I noticed that it does not actually "reverse" anything, but is it necessary?
Thanks, again. 
-----Original Message-----
From: David L Carlson <[hidden email]>
To: [hidden email] <[hidden email]>; [hidden email] <[hidden email]>
Sent: Fri, Apr 5, 2019 12:59 pm
Subject: RE: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

The two lines can have different x coordinates. We are defining the boundaries of the polygon by tracing clockwise or counterclockwise around it, e.g. In this example we follow the top line from beginning to end and then follow the bottom line from end back to beginning. The function automatically connects the beginning and end points to define the polygon and then shades it. The x values do not need to be the same for the two lines, e.g.

lines(x1, y1)
lines(x2, y2)
polygon(c(x1, rev(x2)), c(y1, rev(y1)))

Things are a bit trickier if you want the polygon to start and end at the edge of the plot window. You can get the coordinates of the plot window boundaries with the following:

par("usr")
# [1]  0.04 25.96 -0.20  5.20

This gives the left lower corner (.04, -.02) and the upper right corner (25.96, 5.20). We only need the x values: .04 on the left and 25.96 on the right. We insert x-coordinates into the polygon, but we have to compute the y-axis values at .04 and 25.96 using the formulae you are using to define those lines.

x <- 1:25
y1 <- sqrt(x)
y2 <- x/6
y3 <- y2 + .5
plot(x, y1, typ="l", ylim=c(0, 5))
lines(x, y2)
lines(x, y3)
greena <- rgb(0, 255/255, 0, .5)

edge <- par("usr")
polygon(c(edge[1], x, 25.96, 25.96, rev(x), edge[1]) ,
    c(sqrt(edge[1]), y1, sqrt(edge[2]), edge[2]/6, rev(y2), edge[1]/6),
    col=greena)

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352


From: [hidden email] <[hidden email]>
Sent: Friday, April 5, 2019 11:26 AM
To: David L Carlson <[hidden email]>; [hidden email]
Subject: Re: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

Hi David,
Okay, I just gave this a try and was able to generate a shaded area using my existing code. But just to make sure that I am using this correctly, I did this using the following

plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
polygon(c(get,rev(get)), c(Hope2,rev(Hope9)),col="green")   #this is what I just tried



The polygon line that I just added there - am I telling R to shade everything from Hope2 to Hope9?

In this case, I have different x and y values for each "lines". Do I need to specify the polygon function each time I use the "lines" function in order to capture "all" of the green curves?

Thanks, again,

-----Original Message-----
From: David L Carlson <mailto:[hidden email]>
To: mailto:[hidden email] <mailto:[hidden email]>; mailto:[hidden email] <mailto:[hidden email]>
Sent: Fri, Apr 5, 2019 11:53 am
Subject: RE: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R
You can use the polygon() function (?polygon). Here's a simple reproducible example:

x <- 1:25
y1 <- sqrt(x)
y2 <- x/6
plot(x, y1, typ="l", ylim=c(0, 5))
lines(x, y2)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col="green")
lines(x, y2+.5)

If there are lines hidden by the shading, plot them after the polygon or add an alpha value to the fill color you are using, e.g.

greena <- rgb(0, 255/255, 0, .5)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col=greena)

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352

-----Original Message-----
From: R-sig-Geo <mailto:[hidden email]> On Behalf Of rain1290--- via R-sig-Geo
Sent: Friday, April 5, 2019 10:22 AM
To: mailto:[hidden email]
Subject: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

Hello there,
I have many line plots on a single graph and would like to somehow create a shaded area between them so that it looks neater and easier to interpret. Something like this:
http://jvoigts.scripts.mit.edu/blog/assets/plot_shaded_pretty.png

I heard that the function "geom_ribbon" was ideal to do this, but is it possible to incorporate this in the basic "plot" command? If so, how?
This is currently what I have done with the basic plotting, which works fine:
plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
So, in this plot, I have 5 green line plots in one. The idea would be shade these in green, from the lower to upper bounds. Is this possible?
Any assistance would be extremely valuable!

    [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
mailto:[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: Creating a color shade between upper and lower bounds of a line plot in R

Fri, 04/05/2019 - 11:59
The two lines can have different x coordinates. We are defining the boundaries of the polygon by tracing clockwise or counterclockwise around it, e.g. In this example we follow the top line from beginning to end and then follow the bottom line from end back to beginning. The function automatically connects the beginning and end points to define the polygon and then shades it. The x values do not need to be the same for the two lines, e.g.

lines(x1, y1)
lines(x2, y2)
polygon(c(x1, rev(x2)), c(y1, rev(y1)))

Things are a bit trickier if you want the polygon to start and end at the edge of the plot window. You can get the coordinates of the plot window boundaries with the following:

par("usr")
# [1]  0.04 25.96 -0.20  5.20

This gives the left lower corner (.04, -.02) and the upper right corner (25.96, 5.20). We only need the x values: .04 on the left and 25.96 on the right. We insert x-coordinates into the polygon, but we have to compute the y-axis values at .04 and 25.96 using the formulae you are using to define those lines.

x <- 1:25
y1 <- sqrt(x)
y2 <- x/6
y3 <- y2 + .5
plot(x, y1, typ="l", ylim=c(0, 5))
lines(x, y2)
lines(x, y3)
greena <- rgb(0, 255/255, 0, .5)

edge <- par("usr")
polygon(c(edge[1], x, 25.96, 25.96, rev(x), edge[1]) ,
     c(sqrt(edge[1]), y1, sqrt(edge[2]), edge[2]/6, rev(y2), edge[1]/6),
     col=greena)

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352


From: [hidden email] <[hidden email]>
Sent: Friday, April 5, 2019 11:26 AM
To: David L Carlson <[hidden email]>; [hidden email]
Subject: Re: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

Hi David,
Okay, I just gave this a try and was able to generate a shaded area using my existing code. But just to make sure that I am using this correctly, I did this using the following

plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
polygon(c(get,rev(get)), c(Hope2,rev(Hope9)),col="green")   #this is what I just tried



The polygon line that I just added there - am I telling R to shade everything from Hope2 to Hope9?

In this case, I have different x and y values for each "lines". Do I need to specify the polygon function each time I use the "lines" function in order to capture "all" of the green curves?

Thanks, again,

-----Original Message-----
From: David L Carlson <mailto:[hidden email]>
To: mailto:[hidden email] <mailto:[hidden email]>; mailto:[hidden email] <mailto:[hidden email]>
Sent: Fri, Apr 5, 2019 11:53 am
Subject: RE: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R
You can use the polygon() function (?polygon). Here's a simple reproducible example:

x <- 1:25
y1 <- sqrt(x)
y2 <- x/6
plot(x, y1, typ="l", ylim=c(0, 5))
lines(x, y2)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col="green")
lines(x, y2+.5)

If there are lines hidden by the shading, plot them after the polygon or add an alpha value to the fill color you are using, e.g.

greena <- rgb(0, 255/255, 0, .5)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col=greena)

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352

-----Original Message-----
From: R-sig-Geo <mailto:[hidden email]> On Behalf Of rain1290--- via R-sig-Geo
Sent: Friday, April 5, 2019 10:22 AM
To: mailto:[hidden email]
Subject: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

Hello there,
I have many line plots on a single graph and would like to somehow create a shaded area between them so that it looks neater and easier to interpret. Something like this:
http://jvoigts.scripts.mit.edu/blog/assets/plot_shaded_pretty.png

I heard that the function "geom_ribbon" was ideal to do this, but is it possible to incorporate this in the basic "plot" command? If so, how?
This is currently what I have done with the basic plotting, which works fine:
plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
So, in this plot, I have 5 green line plots in one. The idea would be shade these in green, from the lower to upper bounds. Is this possible?
Any assistance would be extremely valuable!

    [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
mailto:[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: Creating a color shade between upper and lower bounds of a line plot in R

Fri, 04/05/2019 - 11:26
Hi David,

Okay, I just gave this a try and was able to generate a shaded area using my existing code. But just to make sure that I am using this correctly, I did this using the following
plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
polygon(c(get,rev(get)), c(Hope2,rev(Hope9)),col="green")   #this is what I just tried


The polygon line that I just added there - am I telling R to shade everything from Hope2 to Hope9?
In this case, I have different x and y values for each "lines". Do I need to specify the polygon function each time I use the "lines" function in order to capture "all" of the green curves?
Thanks, again,
-----Original Message-----
From: David L Carlson <[hidden email]>
To: [hidden email] <[hidden email]>; [hidden email] <[hidden email]>
Sent: Fri, Apr 5, 2019 11:53 am
Subject: RE: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

You can use the polygon() function (?polygon). Here's a simple reproducible example:

x <- 1:25
y1 <- sqrt(x)
y2 <- x/6
plot(x, y1, typ="l", ylim=c(0, 5))
lines(x, y2)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col="green")
lines(x, y2+.5)

If there are lines hidden by the shading, plot them after the polygon or add an alpha value to the fill color you are using, e.g.

greena <- rgb(0, 255/255, 0, .5)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col=greena)

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352

-----Original Message-----
From: R-sig-Geo <[hidden email]> On Behalf Of rain1290--- via R-sig-Geo
Sent: Friday, April 5, 2019 10:22 AM
To: [hidden email]
Subject: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

Hello there,
I have many line plots on a single graph and would like to somehow create a shaded area between them so that it looks neater and easier to interpret. Something like this:
http://jvoigts.scripts.mit.edu/blog/assets/plot_shaded_pretty.png

I heard that the function "geom_ribbon" was ideal to do this, but is it possible to incorporate this in the basic "plot" command? If so, how?
This is currently what I have done with the basic plotting, which works fine:
plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
So, in this plot, I have 5 green line plots in one. The idea would be shade these in green, from the lower to upper bounds. Is this possible?
Any assistance would be extremely valuable!

    [[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: Creating a color shade between upper and lower bounds of a line plot in R

Fri, 04/05/2019 - 10:53
You can use the polygon() function (?polygon). Here's a simple reproducible example:

x <- 1:25
y1 <- sqrt(x)
y2 <- x/6
plot(x, y1, typ="l", ylim=c(0, 5))
lines(x, y2)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col="green")
lines(x, y2+.5)

If there are lines hidden by the shading, plot them after the polygon or add an alpha value to the fill color you are using, e.g.

greena <- rgb(0, 255/255, 0, .5)
polygon(c(x, rev(x)) , c(y1, rev(y2)), col=greena)

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352

-----Original Message-----
From: R-sig-Geo <[hidden email]> On Behalf Of rain1290--- via R-sig-Geo
Sent: Friday, April 5, 2019 10:22 AM
To: [hidden email]
Subject: [R-sig-Geo] Creating a color shade between upper and lower bounds of a line plot in R

Hello there,
I have many line plots on a single graph and would like to somehow create a shaded area between them so that it looks neater and easier to interpret. Something like this:
http://jvoigts.scripts.mit.edu/blog/assets/plot_shaded_pretty.png

I heard that the function "geom_ribbon" was ideal to do this, but is it possible to incorporate this in the basic "plot" command? If so, how?
This is currently what I have done with the basic plotting, which works fine:
plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
So, in this plot, I have 5 green line plots in one. The idea would be shade these in green, from the lower to upper bounds. Is this possible?
Any assistance would be extremely valuable!

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

Creating a color shade between upper and lower bounds of a line plot in R

Fri, 04/05/2019 - 10:21
Hello there,
I have many line plots on a single graph and would like to somehow create a shaded area between them so that it looks neater and easier to interpret. Something like this:
http://jvoigts.scripts.mit.edu/blog/assets/plot_shaded_pretty.png

I heard that the function "geom_ribbon" was ideal to do this, but is it possible to incorporate this in the basic "plot" command? If so, how?
This is currently what I have done with the basic plotting, which works fine:
plot(get, Hope2, type = "l",col = "green", lwd = "3", xlab="Cumulative CO2 emissions (TtC)", ylab = "One-day maximum precipitation (mm/day)", main = "One-day maximum precipitation for Random Area for CanESM2 scenarios")lines(IPSL, Hope6, type = "l", lwd = "3", col = "green")
lines(MIROC, Hope7, type = "l", lwd = "3", col = "green")
lines(subsetprime, Hope8, type = "l", lwd = "3", col = "green")
lines(MPI, Hope9, type =" l", lwd = "3", col = "green")
So, in this plot, I have 5 green line plots in one. The idea would be shade these in green, from the lower to upper bounds. Is this possible?
Any assistance would be extremely valuable!

        [[alternative HTML version deleted]]

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

Re: Common Factors test in R

Fri, 04/05/2019 - 08:11
Hello,
I would like to perform a common factor test, conducted using likelihood ratio test on spatial error model and spatial durbin model (both in panel format).
I have not found a command in R, that would help me to conduct the test. Hence, I am trying to perform the test manuallly using the splm package and data available in the splm and plm packages.
Here is my code:
library(splm)
library(plm)

data(Produc, package="plm")
data(usaww)

Produc <- pdata.frame(Produc, index = c("state", "year"))

usaww<- mat2listw(usaww, style="W")

Produc$slagUnemp <- slag(Produc$unemp, listw = usaww)

durbin <- spml(gsp~unemp + slagUnemp,
                data=Produc, listw=usaww,  effect = "twoways", model="within", lag=TRUE, spatial.error = "none", quiet = FALSE)
spem  <- spml(gsp ~ unemp,
                data=Produc, listw=usaww,  effect = "twoways", model="within", lag=FALSE, spatial.error = "b", quiet = FALSE)
Is it correct to take the last value of the function from the console output to compute the likelihood ratio?
Meaning, in this example, to calculate 2*(function(non-nested) - function(nested)) = 2*(10261.79 - 10255.74) = 12.1?
If it is correct, how can I compute the significance values?
If it is not correct, is there a better way to compute this?

Autorom tejto spr�vy elektronickej po�ty je Martin Hul�nyi. T�to spr�va je ur�en� v�lu�ne jej adres�tovi. Inform�cie a �daje, ktor� s� v nej uveden�, alebo ktor� s� obsiahnut� v jej prilo�en�ch s�boroch, m��u by� inform�ciami alebo �dajmi chr�nen�mi pod�a platn�ch pr�vnych predpisov v Slovenskej republike. V pr�pade, ak nie ste ur�en� ako prij�mate� tejto spr�vy alebo jeho opr�vnen� z�stupca, upozor�ujeme V�s, �e inform�cie a �daje v nej uveden� nie ste opr�vnen� sprac�va�, ani ich spr�stupni� alebo poskytn�� tretej osobe alebo ich zverejni�. Ak ste nedopatren�m prijali alebo zachytili tuto spr�vu elektronickej po�ty, dovo�ujeme si V�s po�iada�, aby ste ju zaslali sp� na elektronick� adresu jej odosielate�a, a aby ste ju n�sledne zmazali zo svojho po��ta�a a z Va�ej schr�nky elektronickej po�ty. Odosielate� tejto spr�vy nenesie zodpovednos� za �kody sp�soben� nespr�vnym pou�it�m tejto spr�vy elektronickej po�ty a jej pr�loh.

The author of this e-mail message is Martin Hul�nyi. This message is intended solely for the recipient. The information and data contained therein or contained in its enclosed files may be information or data protected under applicable laws in the Slovak Republic. If you are not designated as the recipient of this message or its authorized representative, we would like to inform you that the information and data contained therein are not authorized to process or make them available to third parties or to disclose them. If you have received or downloaded this e-mail message accidentally, we ask you to send it back to the e-mail address of the sender and then delete it from your computer and from your e-mail. The sender of this message is not responsible for damages caused by incorrect use of this e-mail message and its attachments.
[eco.jpg]       Pred vytla�en�m tohto mailu zv�te pros�m vplyv na �ivotn� prostredie. �akujeme.
Please consider the environment before printing this e-mail. Thanks

        [[alternative HTML version deleted]]


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

Mosaic

Thu, 04/04/2019 - 09:45
Hi,
I am working with time series Landsat data and creating mosaic images.
What would you offer me to create seamless, color balanced, and smooth mosaic images from multiple Landsat 8 images?

Thanks,

--
Fatih Kara

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

Re: Match coordinates with NUTS3 regions

Thu, 04/04/2019 - 00:53
Dear Michele,
be sure that the order of the coordinates is  the same for the two objects.
proj4string(EU_NUTS29) is starting with "+proj=longlat" (long, then
lat), while coordinates(awards_with_coordinates_first2500) is "~g_lat +
g_lon" (lat, then lon).
HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences

2019.04.04. 1:47 keltezéssel, Michele Molteni írta:
> Hi everyone,
>
> I have a .dta file with coordinates (WGS 84) and I would like to match them
> to NUTS3 regions. I followed indications from a previous post but I get no
> match even though the coordinates fall within NUTS3 bboxes. Could you
> please be so kind to help me understand what I am doing wrong? This is what
> I have so far:
>
>> library(rgdal)
>>
>> ## Download Administrative Level data from EuroStat
>>
>> temp <- tempfile(fileext = ".zip")
>> download.file("
> https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-60m.shp.zip",
> temp)
> trying URL '
> https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-60m.shp.zip
> '
> Content type 'application/zip' length 2330060 bytes (2.2 MB)
> downloaded 2.2 MB
>> unzip(temp)
>>
>> ##Read data after having unzipped it
>> EU_NUTS <- readOGR(dsn = "C:/Users/User/Documents/TED_project/Replicating
> paper/Shapefile", layer="NUTS_RG_60M_2013_4326_LEVL_3")
> OGR data source with driver: ESRI Shapefile
> Source: "C:\Users\User\Documents\TED_project\Replicating paper\Shapefile",
> layer: "NUTS_RG_60M_2013_4326_LEVL_3"
> with 1480 features
> It has 5 fields
>> ##Consider only EU29 countries (1361 NUTS3)
>> EU_NUTS29 <- EU_NUTS[which(EU_NUTS@data$CNTR_CODE %in% c('BE', 'BG',
> 'CZ', 'DK', 'DE', 'EE', 'IE', 'EL', 'ES', 'FR', 'HR', 'IT', 'CY', 'LV',
> 'LT', 'LU', 'HU', 'MT', 'NL', 'AT', 'PL', 'PT', 'RO', 'SI', 'SK', 'FI',
> 'SE', 'UK', 'NO')), ]
>> ##Import data contaning lat and lon and consider first 2500 observations
>> awards_with_coordinates_first2500 <- read_dta("~/TED_project/Replicating
> paper/awards_with_coordinates_first2500.dta")
>> awards_with_coordinates_first2500 <- awards_with_coordinates_first2
> 500[(1:2500),]
>> ## Here is the projection for map_nuts data
>> proj4string(EU_NUTS29)
> [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
>> ##Transform g_lat and g_lon in numeric
>> awards_with_coordinates_first2500$g_lat <- as.numeric(as.character(awards
> _with_coordinates_first2500$g_lat))
>> awards_with_coordinates_first2500$g_lon <- as.numeric(as.character(awards
> _with_coordinates_first2500$g_lon))
>> ## Make a spatial points dataframe from point_data
>> ## Assuming the projection to be longlat on the WGS84 datum
>> ## First define the coordinates
>> coordinates(awards_with_coordinates_first2500) <- ~g_lat + g_lon
>>
>> ## Then the projection
>> proj4string(awards_with_coordinates_first2500) <- CRS("+proj=longlat
> +datum=WGS84")
>> ## Now transform the projection of awards_with_coordinates_first2500 to
> that of EU_NUTS29
>> awards_with_coordinates_first2500 <- spTransform(awards_with_coordinates_first2500,
> proj4string(EU_NUTS29))
>> ## Use the over() function to get the regions for the points
>> over(awards_with_coordinates_first2500,EU_NUTS29,by.id=)
>      NUTS_ID LEVL_CODE CNTR_CODE NUTS_NAME  FID
> 1      <NA>        NA      <NA>      <NA> <NA>
> 2      <NA>        NA      <NA>      <NA> <NA>
>
> 200    <NA>        NA      <NA>      <NA> <NA>
>   [ reached 'max' / getOption("max.print") -- omitted 2300 rows ]
>> bbox(awards_with_coordinates_first2500)
>              min      max
> g_lat 34.794891 51.46780
> g_lon -3.018564 33.63785
>> bbox(EU_NUTS29)
>        min    max
> x -61.806 55.826
> y -21.350 71.127
>> ##The points seem to be all in Africa
>> plot(EU_NUTS29)
>> plot(awards_with_coordinates_first2500, add = TRUE, col = "red", pch = 19)
> Thanks in advance.
>
> Cheers,
> Michele.
>
> [[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

Match coordinates with NUTS3 regions

Wed, 04/03/2019 - 18:47
Hi everyone,

I have a .dta file with coordinates (WGS 84) and I would like to match them
to NUTS3 regions. I followed indications from a previous post but I get no
match even though the coordinates fall within NUTS3 bboxes. Could you
please be so kind to help me understand what I am doing wrong? This is what
I have so far:

> library(rgdal)
>
> ## Download Administrative Level data from EuroStat
>
> temp <- tempfile(fileext = ".zip")
> download.file("
https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-60m.shp.zip",
temp)
trying URL '
https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-60m.shp.zip
'
Content type 'application/zip' length 2330060 bytes (2.2 MB)
downloaded 2.2 MB
>
> unzip(temp)
>
> ##Read data after having unzipped it
> EU_NUTS <- readOGR(dsn = "C:/Users/User/Documents/TED_project/Replicating
paper/Shapefile", layer="NUTS_RG_60M_2013_4326_LEVL_3")
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\User\Documents\TED_project\Replicating paper\Shapefile",
layer: "NUTS_RG_60M_2013_4326_LEVL_3"
with 1480 features
It has 5 fields
>
> ##Consider only EU29 countries (1361 NUTS3)
> EU_NUTS29 <- EU_NUTS[which(EU_NUTS@data$CNTR_CODE %in% c('BE', 'BG',
'CZ', 'DK', 'DE', 'EE', 'IE', 'EL', 'ES', 'FR', 'HR', 'IT', 'CY', 'LV',
'LT', 'LU', 'HU', 'MT', 'NL', 'AT', 'PL', 'PT', 'RO', 'SI', 'SK', 'FI',
'SE', 'UK', 'NO')), ]
>
> ##Import data contaning lat and lon and consider first 2500 observations
> awards_with_coordinates_first2500 <- read_dta("~/TED_project/Replicating
paper/awards_with_coordinates_first2500.dta")
> awards_with_coordinates_first2500 <- awards_with_coordinates_first2
500[(1:2500),]
>
> ## Here is the projection for map_nuts data
> proj4string(EU_NUTS29)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
>
> ##Transform g_lat and g_lon in numeric
> awards_with_coordinates_first2500$g_lat <- as.numeric(as.character(awards
_with_coordinates_first2500$g_lat))
> awards_with_coordinates_first2500$g_lon <- as.numeric(as.character(awards
_with_coordinates_first2500$g_lon))
>
> ## Make a spatial points dataframe from point_data
> ## Assuming the projection to be longlat on the WGS84 datum
> ## First define the coordinates
> coordinates(awards_with_coordinates_first2500) <- ~g_lat + g_lon
>
> ## Then the projection
> proj4string(awards_with_coordinates_first2500) <- CRS("+proj=longlat
+datum=WGS84")
>
> ## Now transform the projection of awards_with_coordinates_first2500 to
that of EU_NUTS29
> awards_with_coordinates_first2500 <- spTransform(awards_with_coordinates_first2500,
proj4string(EU_NUTS29))
>
> ## Use the over() function to get the regions for the points
> over(awards_with_coordinates_first2500,EU_NUTS29,by.id=)
    NUTS_ID LEVL_CODE CNTR_CODE NUTS_NAME  FID
1      <NA>        NA      <NA>      <NA> <NA>
2      <NA>        NA      <NA>      <NA> <NA>

200    <NA>        NA      <NA>      <NA> <NA>
 [ reached 'max' / getOption("max.print") -- omitted 2300 rows ]
>
> bbox(awards_with_coordinates_first2500)
            min      max
g_lat 34.794891 51.46780
g_lon -3.018564 33.63785
> bbox(EU_NUTS29)
      min    max
x -61.806 55.826
y -21.350 71.127
> ##The points seem to be all in Africa
> plot(EU_NUTS29)
> plot(awards_with_coordinates_first2500, add = TRUE, col = "red", pch = 19)

Thanks in advance.

Cheers,
Michele.

        [[alternative HTML version deleted]]

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

Creating shades around a cluster of line plots

Wed, 04/03/2019 - 11:18
Hi there,
I am currently working with a line plot that would contain a large number of separate plots (which does make it difficult to read). What I would like to do instead is to somehow create a light-colored shade that captures the range of individual line plots. Would that be possible? This is what I have for my plot:
plot(get,Hope2, type="l",col="red", lwd="3", xlab="Cumulative CO2 emissions (TtC)", ylab="One-day maximum precipitation (mm/day)", main="One-day maximum precipitation for South Sudan for CanESM2 scenarios")
lines(get2.teratons, Hope3, type="l", lwd="3", col="red")
lines(get3.teratons, Hope4, type="l", lwd="3", col="red")
lines(get4.teratons, Hope5, type="l", lwd="3", col="red")
So, this gives 4 separate red lines on the same plot (I am likely going to place up to 10 lines, so you can imagine how cluttered that would look without shading them in the background!). Now, let's say that I wanted to create a light red shade that runs from the upper to lower red curves. How should I go about doing that? The first line corresponds to the upper line, and the last line is the lower line.
Thanks,

        [[alternative HTML version deleted]]

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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 19:49
SUCCESS! Thank you so, so much, Stephen! Indeed, that worked out well, and dim(subset) shows:
64   128    90
I am very, very grateful for this! Thanks to the others, too, for their assistance with this!

-----Original Message-----
From: Stephen Stewart <[hidden email]>
To: Rolf Turner <[hidden email]>
Cc: rain1290 <[hidden email]>; R-sig-geo Mailing List <[hidden email]>
Sent: Sat, Mar 30, 2019 8:16 pm
Subject: Re: [R-sig-Geo] [FORGED] Modifying the length of a matrix variable

The last few messages provided important information which was previously absent. Try:

subset = Model4[[1:90]]

This will subset the brick to layers 1 through 90.

I would also suggest some further reading around the raster package and NetCDF files (e.g. the ncdf4 package) would be useful to you.

On Sun, 31 Mar 2019 at 11:06, Rolf Turner <[hidden email]> wrote:


On 31/03/19 12:56 PM, [hidden email] wrote:

> Model4 <- brick("MaxPrecCCCMACanESM2rcp45.nc", var="onedaymax")
>
>
> That is how Model4 is derived.
>
> When trying class(Model4), I receive:
>
> [1] "RasterBrick" attr(,"package") [1] "raster"
>
> **//___^
>
> Meanwhile, I will check on Google to see what I come up with in terms of
> your suggestion. :)
I know nothing about rasters, the brick() function, or the raster
package, so include me out at this stage.

Others on the list may be able to help you.  Particular if you can force
yourself to ask a *focussed* question.

cheers,

Rolf

--
Honorary Research Fellow
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


        [[alternative HTML version deleted]]

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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 18:16
The last few messages provided important information which was previously
absent. Try:

subset = Model4[[1:90]]

This will subset the brick to layers 1 through 90.

I would also suggest some further reading around the raster package and
NetCDF files (e.g. the ncdf4 package) would be useful to you.

On Sun, 31 Mar 2019 at 11:06, Rolf Turner <[hidden email]> wrote:

>
> On 31/03/19 12:56 PM, [hidden email] wrote:
>
> > Model4 <- brick("MaxPrecCCCMACanESM2rcp45.nc", var="onedaymax")
> >
> >
> > That is how Model4 is derived.
> >
> > When trying class(Model4), I receive:
> >
> > [1] "RasterBrick" attr(,"package") [1] "raster"
> >
> > **//___^
> >
> > Meanwhile, I will check on Google to see what I come up with in terms of
> > your suggestion. :)
>
> I know nothing about rasters, the brick() function, or the raster
> package, so include me out at this stage.
>
> Others on the list may be able to help you.  Particular if you can force
> yourself to ask a *focussed* question.
>
> cheers,
>
> Rolf
>
> --
> Honorary Research Fellow
> 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
>
        [[alternative HTML version deleted]]

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

PhD Student
School of Ecosystem and Forest Sciences
Faculty of Science
University of Melbourne

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 18:06
Hi,

It looks like you are accessing a NetCDF file or it's kin.  Like Rolf, I'm challenged to understand from your description what you are after.  Assuming you do indeed have a NetCDF file and are using the ncdf4 package, you can retrieve the first 90 layers as shown below using the varid, start and count arguments to ncvar_get().  Obviously, this is untested since your data isn't available to us.

library(ncdf4)
filename = "/path/to/foo.nc"
ncObject <- nc_open(filename)
m <- ncvar_get(ncObject varid = "onedaymax", start = c(1,1,1), count = c(-1,-1, 90))
nc_close(ncObject)


See more details and examples here... https://www.rdocumentation.org/packages/ncdf4/versions/1.16.1/topics/ncvar_get


You might also consider using the raster package which will retrieve multilayer georeferenced rasters for you from the same NetCDF file.

library(raster)
B <- brick(filename, varname = "onedaymax")
B <- B[[-c(91:95)]]    # <- drop the last layers 91, 92, ...95

See the docs here... https://www.rdocumentation.org/packages/raster/versions/2.8-19/topics/brick

Cheers,
Ben

> On Mar 30, 2019, at 7:30 PM, rain1290--- via R-sig-Geo <[hidden email]> wrote:
>
> Yes, I reproduced the example above and it works just fine (and is what I want!!), but I cannot see why it does not work with my data, as it is a 3-dimensional array (latitude, longitude and time).
> This is what comes from print(Model4):
> 3 variables (excluding dimension variables):
>        double onedaymax[lon,lat,time]   (Contiguous storage)  
>            units: mm/day
>        double fivedaymax[lon,lat,time]   (Contiguous storage)  
>            units: mm/day
>        short Year[time]   (Contiguous storage)  
>
>     3 dimensions:
>        time  Size:95
>        lat  Size:64
>            units: degree North
>        lon  Size:128
>            units: degree East
> I reviewed it over and over again, but I cannot see why this would not work?
> Thanks,
>
>
> -----Original Message-----
> From: Rolf Turner <[hidden email]>
> To: rain1290 <[hidden email]>
> Cc: r-sig-geo <[hidden email]>
> Sent: Sat, Mar 30, 2019 6:49 pm
> Subject: Re: [FORGED] [R-sig-Geo] Modifying the length of a matrix variable
>
>
> On 31/03/19 11:14 AM, [hidden email] wrote:
>
>> Hi Rolf (and others),
>>
>> I tried your suggestion, but when I used dim(Model4.chopped), it still
>> shows 95 layers, as shown below:
>>
>> 8192    95
>>
>> I also find that the total number of cells is rather low for that many
>> layers. I started with 778240 cells over 95 layers.
>
> Well then you're doing something wrong, or there is something that you
> haven't told us.
>
> E.g.:
>> junk <- array(runif(64*128*95),dim=c(64,128,95))
>> junk.chopped <- junk[,,1:90]
>> dim(junk)
>> [1]  64 128  95
>> dim(junk.chopped)
>> [1]  64 128  90
>
> Perhaps Model.4 has some structure other than that of an array.
> (Originally you said it was a matrix.)
>
> You really need to get your terminology and ideas *clear* in order to
> have any hope of receiving useful advice.
>
> I have no idea what you are on about in respect of "the number of
> cells".  My mind-reading machine is in the repair shop.  I strongly
> suspect that your thoughts are confused.
>
> cheers,
>
> Rolf
>
> --
> Honorary Research Fellow
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

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

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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 18:01

On 31/03/19 12:56 PM, [hidden email] wrote:

> Model4 <- brick("MaxPrecCCCMACanESM2rcp45.nc", var="onedaymax")
>
>
> That is how Model4 is derived.
>
> When trying class(Model4), I receive:
>
> [1] "RasterBrick" attr(,"package") [1] "raster"
>
> **//___^
>
> Meanwhile, I will check on Google to see what I come up with in terms of
> your suggestion. :)
I know nothing about rasters, the brick() function, or the raster
package, so include me out at this stage.

Others on the list may be able to help you.  Particular if you can force
yourself to ask a *focussed* question.

cheers,

Rolf

--
Honorary Research Fellow
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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 17:56
Model4 <- brick("MaxPrecCCCMACanESM2rcp45.nc", var="onedaymax")
That is how Model4 is derived. 
When trying class(Model4), I receive:
[1] "RasterBrick"
attr(,"package")
[1] "raster"

Meanwhile, I will check on Google to see what I come up with in terms of your suggestion. :)
-----Original Message-----
From: Rolf Turner <[hidden email]>
To: rain1290 <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Sat, Mar 30, 2019 7:43 pm
Subject: Re: [FORGED] [R-sig-Geo] Modifying the length of a matrix variable


On 31/03/19 12:30 PM, [hidden email] wrote:

> Yes, I reproduced the example above and it works just fine (and is what
> I want!!), but I cannot see why it does not work with my data, as it is
> a 3-dimensional array (latitude, longitude and time).
>
> This is what comes from print(Model4):
>
> 3 variables (excluding dimension variables): double
> onedaymax[lon,lat,time] (Contiguous storage) units: mm/day double
> fivedaymax[lon,lat,time] (Contiguous storage) units: mm/day short
> Year[time] (Contiguous storage) 3 dimensions: time Size:95 lat Size:64
> units: degree North lon Size:128 units: degree East
>
>
> I reviewed it over and over again, but I cannot see why this would not work?
Psigh!  Clearly Model4 is *not* an array!!!  It is an object of some
"specialised" class (for which there is specialised print() method).  I
have no idea what that class might be, but *you can tell.  What does

    class(Model4)

return?

Where did this "Model4" object come from?  What are you trying to *do*?

You might be able to get somewhere by searching (e.g. via Google) on
"subsetting objects of class melvin" where "melvin" is what is returned
by "class(Model4)".

Doing

    str(Model4)

could be enlightening (but given your stubborn refusal to acquire
insight into the workings of R, I am not optimistic).

This is not magic or religion.  You need to *understand* what you are
dealing with, and proceed rationally.  Don't just hammer and hope.

cheers,

Rolf

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

        [[alternative HTML version deleted]]

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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 17:43

On 31/03/19 12:30 PM, [hidden email] wrote:

> Yes, I reproduced the example above and it works just fine (and is what
> I want!!), but I cannot see why it does not work with my data, as it is
> a 3-dimensional array (latitude, longitude and time).
>
> This is what comes from print(Model4):
>
> 3 variables (excluding dimension variables): double
> onedaymax[lon,lat,time] (Contiguous storage) units: mm/day double
> fivedaymax[lon,lat,time] (Contiguous storage) units: mm/day short
> Year[time] (Contiguous storage) 3 dimensions: time Size:95 lat Size:64
> units: degree North lon Size:128 units: degree East
>
>
> I reviewed it over and over again, but I cannot see why this would not work?
Psigh!  Clearly Model4 is *not* an array!!!  It is an object of some
"specialised" class (for which there is specialised print() method).  I
have no idea what that class might be, but *you can tell.  What does

    class(Model4)

return?

Where did this "Model4" object come from?  What are you trying to *do*?

You might be able to get somewhere by searching (e.g. via Google) on
"subsetting objects of class melvin" where "melvin" is what is returned
by "class(Model4)".

Doing

     str(Model4)

could be enlightening (but given your stubborn refusal to acquire
insight into the workings of R, I am not optimistic).

This is not magic or religion.  You need to *understand* what you are
dealing with, and proceed rationally.  Don't just hammer and hope.

cheers,

Rolf

--
Honorary Research Fellow
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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 17:30
Yes, I reproduced the example above and it works just fine (and is what I want!!), but I cannot see why it does not work with my data, as it is a 3-dimensional array (latitude, longitude and time). 
This is what comes from print(Model4):
3 variables (excluding dimension variables):
        double onedaymax[lon,lat,time]   (Contiguous storage)  
            units: mm/day
        double fivedaymax[lon,lat,time]   (Contiguous storage)  
            units: mm/day
        short Year[time]   (Contiguous storage)  

     3 dimensions:
        time  Size:95
        lat  Size:64
            units: degree North
        lon  Size:128
            units: degree East
I reviewed it over and over again, but I cannot see why this would not work?
Thanks, 


-----Original Message-----
From: Rolf Turner <[hidden email]>
To: rain1290 <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Sat, Mar 30, 2019 6:49 pm
Subject: Re: [FORGED] [R-sig-Geo] Modifying the length of a matrix variable


On 31/03/19 11:14 AM, [hidden email] wrote:

> Hi Rolf (and others),
>
> I tried your suggestion, but when I used dim(Model4.chopped), it still
> shows 95 layers, as shown below:
>
> 8192    95
>
> I also find that the total number of cells is rather low for that many
> layers. I started with 778240 cells over 95 layers.

Well then you're doing something wrong, or there is something that you
haven't told us.

E.g.:
> junk <- array(runif(64*128*95),dim=c(64,128,95))
> junk.chopped <- junk[,,1:90]
> dim(junk)
> [1]  64 128  95
> dim(junk.chopped)
> [1]  64 128  90

Perhaps Model.4 has some structure other than that of an array.
(Originally you said it was a matrix.)

You really need to get your terminology and ideas *clear* in order to
have any hope of receiving useful advice.

I have no idea what you are on about in respect of "the number of
cells".  My mind-reading machine is in the repair shop.  I strongly
suspect that your thoughts are confused.

cheers,

Rolf

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

        [[alternative HTML version deleted]]

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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 16:49

On 31/03/19 11:14 AM, [hidden email] wrote:

> Hi Rolf (and others),
>
> I tried your suggestion, but when I used dim(Model4.chopped), it still
> shows 95 layers, as shown below:
>
> 8192    95
>
> I also find that the total number of cells is rather low for that many
> layers. I started with 778240 cells over 95 layers.

Well then you're doing something wrong, or there is something that you
haven't told us.

E.g.:
> junk <- array(runif(64*128*95),dim=c(64,128,95))
> junk.chopped <- junk[,,1:90]
> dim(junk)
> [1]  64 128  95
> dim(junk.chopped)
> [1]  64 128  90

Perhaps Model.4 has some structure other than that of an array.
(Originally you said it was a matrix.)

You really need to get your terminology and ideas *clear* in order to
have any hope of receiving useful advice.

I have no idea what you are on about in respect of "the number of
cells".  My mind-reading machine is in the repair shop.  I strongly
suspect that your thoughts are confused.

cheers,

Rolf

--
Honorary Research Fellow
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

Re: [FORGED] Modifying the length of a matrix variable

Sat, 03/30/2019 - 16:14
Hi Rolf (and others),

I tried your suggestion, but when I used dim(Model4.chopped), it still shows 95 layers, as shown below:
8192    95 
I also find that the total number of cells is rather low for that many layers. I started with 778240 cells over 95 layers. 
-----Original Message-----
From: Rolf Turner <[hidden email]>
To: rain1290 <[hidden email]>
Cc: [hidden email] <[hidden email]>
Sent: Sat, Mar 30, 2019 5:57 pm
Subject: Re: [FORGED] [R-sig-Geo] Modifying the length of a matrix variable


On 31/03/19 10:22 AM, [hidden email] wrote:

> Hi Rolf,
>
> My apologies - I meant "layers" as opposed to "length". The goal is to
> reduce the number of layers to 90 (from 95).
>
> dim (Model4) yields:
>
> 64   128   95
>
>
> You can see the 95 there. That is what I would like to reduce to 90, or
> isolate layer 1 to layer 90.
Please keep the list in the set of recipients.  I am not your private
consultant, and furthermore others on the list may be able to provide
better advice than I.  I have CC-ed this message to the list.

To keep only "layers" 1 through 90 you could do:

    Model4.chopped <- Model4[,,1:90]

As I said before, it really is time that you learned something about R
(e.g. by studying a tutorial).

cheers,

Rolf

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

        [[alternative HTML version deleted]]

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

Pages