This is an

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

### Re: [FORGED] 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

_______________________________________________

R-sig-Geo mailing list

[hidden email]

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