Isolating only land areas on a global map for computing averages

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Isolating only land areas on a global map for computing averages

R-sig-geo mailing list
Hi there,
I am interested in calculating precipitation medians globally. However, I only want to isolate land areas to compute the median. I already first created a raster stack, called "RCP1pctCO2median", which contains the median values. There are 138 layers, with each layer representing one year.  This raster stack has the following attributes:
class       : RasterStack
dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
resolution  : 2.8125, 2.789327  (x, y)
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730, 0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687, 0.51857087, 0.41005131, 0.45956529, 0.47497867, ...
max values  :   96.30350,  104.08584,   88.92751,   97.49373,   89.57201,   90.58570,   86.67651,   88.33519,   96.94720,  101.58247,   96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...  

Previously, I was isolating a specific region by specifying a range of longitudes and latitudes to obtain the medians for that region, like this:
expansion1<-expand.grid(103:120, 3:15)lonlataaa<-extract(RCP1pctCO2Median, expansion1)Columnaaa<-colMeans(lonlataaa)

However, with this approach, too much water can mix with land areas, and if I narrow the latitude/longitude range on land, I might miss too much land to compute the median meaningfully.
Therefore, with this data, would it be possible to use an IF/ELSE statement to tell R that if the "center point" of each grid cell happens to fall on land, then it would be considered as land (i.e. that would be TRUE - if not, then FALSE)? Even if a grid cell happens to have water mixed with land, but the center point of the grid is on land, that would be considered land. But can R even tell what is land or water in this case?
Thank you, and I would greatly appreciate any assistance!

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Isolating only land areas on a global map for computing averages

Tom Philippi
The easiest approach would be to create a separate aligned raster layer for
land vs water.  There are plenty of coastline polygons available out there
(e.g., maptools, rworldmap, rnaturalearth packages): choose one in your
raster CRS (or choose one and spTransform() it).  Then, use a grid version
of your raster to extract values from that land/water SpatialPolygons
object.

1: Your idea of extracting the land/water value at each grid cell centroid
gives one estimate.  Instead of TRUE/FALSE, think of the numeric
equivalents 1,0,  then using those as weights for averaging across your
grid cells.
2: A "better" estimate would be to compute the fraction of each grid cell
that is land, then use those fractional [0, 1] values as weights to compute
a weighted average of precipitation over land.  At 2.8deg grid cells, a lot
of heavy rainfall coastal areas would have the grid cell centroid offshore
and be omitted by approach #1.
3: I recommend that you think hard about averaging across cells in Lat Lon
to estimate average precipitation over land.  The actual area of a ~2.8 by
2.8 deg grid cell at the equator is much larger than the same at 70 deg N.
I would spend the extra hour computing the actual area (in km^2) of land in
each of your 8192 grid cells, then using those in a raster as weights for
whatever calculations you do on the raster stack.  [Or you can cheat, as
the area of a grid cell in degrees is a function of only the latitudes, and
your required weights are multiplicative.]

Your mileage may vary...

Tom

On Wed, Nov 6, 2019 at 6:18 PM rain1290--- via R-sig-Geo <
[hidden email]> wrote:

> Hi there,
> I am interested in calculating precipitation medians globally. However, I
> only want to isolate land areas to compute the median. I already first
> created a raster stack, called "RCP1pctCO2median", which contains the
> median values. There are 138 layers, with each layer representing one
> year.  This raster stack has the following attributes:
> class       : RasterStack
> dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
> resolution  : 2.8125, 2.789327  (x, y)
> extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> names       :    layer.1,    layer.2,    layer.3,    layer.4,
> layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,
> layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
> min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730,
> 0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687,
> 0.51857087, 0.41005131, 0.45956529, 0.47497867, ...
> max values  :   96.30350,  104.08584,   88.92751,   97.49373,
> 89.57201,   90.58570,   86.67651,   88.33519,   96.94720,  101.58247,
> 96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...
>
> Previously, I was isolating a specific region by specifying a range of
> longitudes and latitudes to obtain the medians for that region, like this:
> expansion1<-expand.grid(103:120, 3:15)lonlataaa<-extract(RCP1pctCO2Median,
> expansion1)Columnaaa<-colMeans(lonlataaa)
>
> However, with this approach, too much water can mix with land areas, and
> if I narrow the latitude/longitude range on land, I might miss too much
> land to compute the median meaningfully.
> Therefore, with this data, would it be possible to use an IF/ELSE
> statement to tell R that if the "center point" of each grid cell happens to
> fall on land, then it would be considered as land (i.e. that would be TRUE
> - if not, then FALSE)? Even if a grid cell happens to have water mixed with
> land, but the center point of the grid is on land, that would be considered
> land. But can R even tell what is land or water in this case?
> Thank you, and I would greatly appreciate any assistance!
>
>         [[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
Reply | Threaded
Open this post in threaded view
|

Re: Isolating only land areas on a global map for computing averages

R-sig-geo mailing list
Hi Tom and others,

Thank you for your response and suggestions! 
Yes, I loaded and used the maptools package previously to create continents on my world map, among other things. I do think that the easiest approach would be to create a raster layer for land, and then water, with the values that I have. However, my precipitation values are globally distributed - every grid cell has a precipitation value for each year (i.e. each grid cell has 138 values/layers/years). So, if I were to create a raster of only land areas, how would I have my grid cells coincide with the land areas only on that raster?
Also, would it be possible to accomplish this with the raster stack that I already have? If so, is there a way to separate all land/water areas this way using the maptools package?
Thanks, again, and I really appreciate your help!
-----Original Message-----
From: Tom Philippi <[hidden email]>
To: rain1290 <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Thu, Nov 7, 2019 12:44 am
Subject: Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages

The easiest approach would be to create a separate aligned raster layer for land vs water.  There are plenty of coastline polygons available out there (e.g., maptools, rworldmap, rnaturalearth packages): choose one in your raster CRS (or choose one and spTransform() it).  Then, use a grid version of your raster to extract values from that land/water SpatialPolygons object.
1: Your idea of extracting the land/water value at each grid cell centroid gives one estimate.  Instead of TRUE/FALSE, think of the numeric equivalents 1,0,  then using those as weights for averaging across your grid cells.2: A "better" estimate would be to compute the fraction of each grid cell that is land, then use those fractional [0, 1] values as weights to compute a weighted average of precipitation over land.  At 2.8deg grid cells, a lot of heavy rainfall coastal areas would have the grid cell centroid offshore and be omitted by approach #1.3: I recommend that you think hard about averaging across cells in Lat Lon to estimate average precipitation over land.  The actual area of a ~2.8 by 2.8 deg grid cell at the equator is much larger than the same at 70 deg N.  I would spend the extra hour computing the actual area (in km^2) of land in each of your 8192 grid cells, then using those in a raster as weights for whatever calculations you do on the raster stack.  [Or you can cheat, as the area of a grid cell in degrees is a function of only the latitudes, and your required weights are multiplicative.]
Your mileage may vary...
Tom
On Wed, Nov 6, 2019 at 6:18 PM rain1290--- via R-sig-Geo <[hidden email]> wrote:

Hi there,
I am interested in calculating precipitation medians globally. However, I only want to isolate land areas to compute the median. I already first created a raster stack, called "RCP1pctCO2median", which contains the median values. There are 138 layers, with each layer representing one year.  This raster stack has the following attributes:
class       : RasterStack
dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
resolution  : 2.8125, 2.789327  (x, y)
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730, 0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687, 0.51857087, 0.41005131, 0.45956529, 0.47497867, ...
max values  :   96.30350,  104.08584,   88.92751,   97.49373,   89.57201,   90.58570,   86.67651,   88.33519,   96.94720,  101.58247,   96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...  

Previously, I was isolating a specific region by specifying a range of longitudes and latitudes to obtain the medians for that region, like this:
expansion1<-expand.grid(103:120, 3:15)lonlataaa<-extract(RCP1pctCO2Median, expansion1)Columnaaa<-colMeans(lonlataaa)

However, with this approach, too much water can mix with land areas, and if I narrow the latitude/longitude range on land, I might miss too much land to compute the median meaningfully.
Therefore, with this data, would it be possible to use an IF/ELSE statement to tell R that if the "center point" of each grid cell happens to fall on land, then it would be considered as land (i.e. that would be TRUE - if not, then FALSE)? Even if a grid cell happens to have water mixed with land, but the center point of the grid is on land, that would be considered land. But can R even tell what is land or water in this case?
Thank you, and I would greatly appreciate any assistance!

        [[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
Reply | Threaded
Open this post in threaded view
|

Re: Isolating only land areas on a global map for computing averages

Tom Philippi
Rain--
Yes, it is possible to do it with your extant raster stack.  In fact,
pretty much all reasonable approaches will do that.  Anything you do will
create a raster layer with values at every one of your 8192 raster cells:
some values will be 0 or NA.

The trivial answer if you only had a small range of latitude, and small
raster cell sizes relative to your polygon, would be the
raster::rasterize() function to generate a raster mask.

But, with a global raster from +90 to -90, any interpretable answer
averaging over area needs to take into account the different area of a grid
cell at 70 deg latitude vs 0 deg.  See:
https://gis.stackexchange.com/questions/29734/how-to-calculate-area-of-1-x-1-degree-cells-in-a-raster
At that point, you should go ahead and account for fractions of grid cells
over land so your averaging would be over land area.

I'm reticent to give you a complete code solution because without a name,
you may well be a student with this as an assignment.  But, my approach
would be:

create a SpatialPolygonsDataFrame object "gridPoly" from any of your raster
layers via raster::rasterToPolygons()
generate an id value for each of those polygons as row & column indices
from the raster, or as the cell number.
get land polygon SpatialPolygons object "landPoly" with polygons of land
only, not ocean.
create a new SpatialPolygons object from gIntersect::gridPoly, landPoly,
byid = c(TRUE, FALSE))
calculate an area of each polygon in that object via
geosphere::areaPolygon()  {rgeos::gArea()  only works for projected CRS}
create your mask/weights raster layer with either all NA or all 0.
either:
     parse the id values to row & column values.
     use the triplets of row, column, and area to replace the corresponding
NA or 0 values in that mask
or:
     if you used cell numbers, just use the cell numbers and area values in
replacement in that mask

create a second weight raster via raster::area() on one of your raster
layers.
raster multiply your polygon-area and your raster::area values to give the
actual weighs to use.

This still is an approximation, but likely +/- 1-2%.

If this is still complete gibberish to you, either I need more coffee or
you need to consult a good reference on working with spatial data in
general.

On Thu, Nov 7, 2019 at 8:25 AM <[hidden email]> wrote:

> Hi Tom and others,
>
> Thank you for your response and suggestions!
>
> Yes, I loaded and used the maptools package previously to create
> continents on my world map, among other things. I do think that the easiest
> approach would be to create a raster layer for land, and then water, with
> the values that I have. However, my precipitation values are globally
> distributed - every grid cell has a precipitation value for each year (i.e.
> each grid cell has 138 values/layers/years). So, if I were to create a
> raster of only land areas, how would I have my grid cells coincide with the
> land areas only on that raster?
>
> Also, would it be possible to accomplish this with the raster stack that I
> already have? If so, is there a way to separate all land/water areas this
> way using the maptools package?
>
> Thanks, again, and I really appreciate your help!
>
> -----Original Message-----
> From: Tom Philippi <[hidden email]>
> To: rain1290 <[hidden email]>
> Cc: r-sig-geo <[hidden email]>
> Sent: Thu, Nov 7, 2019 12:44 am
> Subject: Re: [R-sig-Geo] Isolating only land areas on a global map for
> computing averages
>
> The easiest approach would be to create a separate aligned raster layer
> for land vs water.  There are plenty of coastline polygons available out
> there (e.g., maptools, rworldmap, rnaturalearth packages): choose one in
> your raster CRS (or choose one and spTransform() it).  Then, use a grid
> version of your raster to extract values from that land/water
> SpatialPolygons object.
>
> 1: Your idea of extracting the land/water value at each grid cell centroid
> gives one estimate.  Instead of TRUE/FALSE, think of the numeric
> equivalents 1,0,  then using those as weights for averaging across your
> grid cells.
> 2: A "better" estimate would be to compute the fraction of each grid cell
> that is land, then use those fractional [0, 1] values as weights to compute
> a weighted average of precipitation over land.  At 2.8deg grid cells, a lot
> of heavy rainfall coastal areas would have the grid cell centroid offshore
> and be omitted by approach #1.
> 3: I recommend that you think hard about averaging across cells in Lat Lon
> to estimate average precipitation over land.  The actual area of a ~2.8 by
> 2.8 deg grid cell at the equator is much larger than the same at 70 deg N.
> I would spend the extra hour computing the actual area (in km^2) of land in
> each of your 8192 grid cells, then using those in a raster as weights for
> whatever calculations you do on the raster stack.  [Or you can cheat, as
> the area of a grid cell in degrees is a function of only the latitudes, and
> your required weights are multiplicative.]
>
> Your mileage may vary...
>
> Tom
>
> On Wed, Nov 6, 2019 at 6:18 PM rain1290--- via R-sig-Geo <
> [hidden email]> wrote:
>
> Hi there,
> I am interested in calculating precipitation medians globally. However, I
> only want to isolate land areas to compute the median. I already first
> created a raster stack, called "RCP1pctCO2median", which contains the
> median values. There are 138 layers, with each layer representing one
> year.  This raster stack has the following attributes:
> class       : RasterStack
> dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
> resolution  : 2.8125, 2.789327  (x, y)
> extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> names       :    layer.1,    layer.2,    layer.3,    layer.4,
> layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,
> layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
> min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730,
> 0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687,
> 0.51857087, 0.41005131, 0.45956529, 0.47497867, ...
> max values  :   96.30350,  104.08584,   88.92751,   97.49373,
> 89.57201,   90.58570,   86.67651,   88.33519,   96.94720,  101.58247,
> 96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...
>
> Previously, I was isolating a specific region by specifying a range of
> longitudes and latitudes to obtain the medians for that region, like this:
> expansion1<-expand.grid(103:120, 3:15)lonlataaa<-extract(RCP1pctCO2Median,
> expansion1)Columnaaa<-colMeans(lonlataaa)
>
> However, with this approach, too much water can mix with land areas, and
> if I narrow the latitude/longitude range on land, I might miss too much
> land to compute the median meaningfully.
> Therefore, with this data, would it be possible to use an IF/ELSE
> statement to tell R that if the "center point" of each grid cell happens to
> fall on land, then it would be considered as land (i.e. that would be TRUE
> - if not, then FALSE)? Even if a grid cell happens to have water mixed with
> land, but the center point of the grid is on land, that would be considered
> land. But can R even tell what is land or water in this case?
> Thank you, and I would greatly appreciate any assistance!
>
>         [[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
Reply | Threaded
Open this post in threaded view
|

Re: Isolating only land areas on a global map for computing averages

R-sig-geo mailing list
Hi Tom and others,

Thank you, once again, for your extensive reply and feedback to this! I really appreciate it!!!
I have some idea as to what you suggested to do, but most of this is still fairly new to me in terms of the procedure. What I am trying to do is actually for personal use, so if you choose to share a complete code demonstrating this process of masking, it could definitely be a useful learning means to better visualize and understand how this is done using my existing raster data. If you choose to, would it be possible to provide some brief comments at each step, so that I can understand what is going on? If not, do you know of any useful references that I can refer to that deals with this procedure? I've been looking, but nothing too specific comes up with what I would like to do, unfortunately.
Also, why would the area of the grid cells at 70-90 degrees latitude differ as compared to the area of lower latitude grid cells if each grid cell is 2.8 x 2.8 degrees? Are you suggesting that I need to derive the cosine of latitude, like the following?
Fncfname <- "MaxPrec5.nc"
FModel <- nc_open(Fncfname)
FPrec  <- brick(Fncfname,var="fivedaymax")
Fsm <- mean(FPrec, na.rm=TRUE)
Fw <- init(FPrec, 'y')
Fw <- cos(Fw*(pi/180))Fx <- Fsm*Fw
Thanks, again, for your extensive help!! 
-----Original Message-----
From: Tom Philippi <[hidden email]>
To: rain1290 <[hidden email]>; tephilippi <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Fri, Nov 8, 2019 11:04 pm
Subject: Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages

Rain--Yes, it is possible to do it with your extant raster stack.  In fact, pretty much all reasonable approaches will do that.  Anything you do will create a raster layer with values at every one of your 8192 raster cells: some values will be 0 or NA.
The trivial answer if you only had a small range of latitude, and small raster cell sizes relative to your polygon, would be the raster::rasterize() function to generate a raster mask.
But, with a global raster from +90 to -90, any interpretable answer averaging over area needs to take into account the different area of a grid cell at 70 deg latitude vs 0 deg.  See: https://gis.stackexchange.com/questions/29734/how-to-calculate-area-of-1-x-1-degree-cells-in-a-rasterAt that point, you should go ahead and account for fractions of grid cells over land so your averaging would be over land area.
I'm reticent to give you a complete code solution because without a name, you may well be a student with this as an assignment.  But, my approach would be:
create a SpatialPolygonsDataFrame object "gridPoly" from any of your raster layers via raster::rasterToPolygons()generate an id value for each of those polygons as row & column indices from the raster, or as the cell number.get land polygon SpatialPolygons object "landPoly" with polygons of land only, not ocean.create a new SpatialPolygons object from gIntersect::gridPoly, landPoly, byid = c(TRUE, FALSE))calculate an area of each polygon in that object via geosphere::areaPolygon()  {rgeos::gArea()  only works for projected CRS}create your mask/weights raster layer with either all NA or all 0.either:      parse the id values to row & column values.     use the triplets of row, column, and area to replace the corresponding NA or 0 values in that mask
or:     if you used cell numbers, just use the cell numbers and area values in replacement in that mask
create a second weight raster via raster::area() on one of your raster layers.
raster multiply your polygon-area and your raster::area values to give the actual weighs to use.
This still is an approximation, but likely +/- 1-2%.
If this is still complete gibberish to you, either I need more coffee or you need to consult a good reference on working with spatial data in general.
On Thu, Nov 7, 2019 at 8:25 AM <[hidden email]> wrote:

Hi Tom and others,

Thank you for your response and suggestions! 
Yes, I loaded and used the maptools package previously to create continents on my world map, among other things. I do think that the easiest approach would be to create a raster layer for land, and then water, with the values that I have. However, my precipitation values are globally distributed - every grid cell has a precipitation value for each year (i.e. each grid cell has 138 values/layers/years). So, if I were to create a raster of only land areas, how would I have my grid cells coincide with the land areas only on that raster?
Also, would it be possible to accomplish this with the raster stack that I already have? If so, is there a way to separate all land/water areas this way using the maptools package?
Thanks, again, and I really appreciate your help!
-----Original Message-----
From: Tom Philippi <[hidden email]>
To: rain1290 <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Thu, Nov 7, 2019 12:44 am
Subject: Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages

The easiest approach would be to create a separate aligned raster layer for land vs water.  There are plenty of coastline polygons available out there (e.g., maptools, rworldmap, rnaturalearth packages): choose one in your raster CRS (or choose one and spTransform() it).  Then, use a grid version of your raster to extract values from that land/water SpatialPolygons object.
1: Your idea of extracting the land/water value at each grid cell centroid gives one estimate.  Instead of TRUE/FALSE, think of the numeric equivalents 1,0,  then using those as weights for averaging across your grid cells.2: A "better" estimate would be to compute the fraction of each grid cell that is land, then use those fractional [0, 1] values as weights to compute a weighted average of precipitation over land.  At 2.8deg grid cells, a lot of heavy rainfall coastal areas would have the grid cell centroid offshore and be omitted by approach #1.3: I recommend that you think hard about averaging across cells in Lat Lon to estimate average precipitation over land.  The actual area of a ~2.8 by 2.8 deg grid cell at the equator is much larger than the same at 70 deg N.  I would spend the extra hour computing the actual area (in km^2) of land in each of your 8192 grid cells, then using those in a raster as weights for whatever calculations you do on the raster stack.  [Or you can cheat, as the area of a grid cell in degrees is a function of only the latitudes, and your required weights are multiplicative.]
Your mileage may vary...
Tom
On Wed, Nov 6, 2019 at 6:18 PM rain1290--- via R-sig-Geo <[hidden email]> wrote:

Hi there,
I am interested in calculating precipitation medians globally. However, I only want to isolate land areas to compute the median. I already first created a raster stack, called "RCP1pctCO2median", which contains the median values. There are 138 layers, with each layer representing one year.  This raster stack has the following attributes:
class       : RasterStack
dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
resolution  : 2.8125, 2.789327  (x, y)
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730, 0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687, 0.51857087, 0.41005131, 0.45956529, 0.47497867, ...
max values  :   96.30350,  104.08584,   88.92751,   97.49373,   89.57201,   90.58570,   86.67651,   88.33519,   96.94720,  101.58247,   96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...  

Previously, I was isolating a specific region by specifying a range of longitudes and latitudes to obtain the medians for that region, like this:
expansion1<-expand.grid(103:120, 3:15)lonlataaa<-extract(RCP1pctCO2Median, expansion1)Columnaaa<-colMeans(lonlataaa)

However, with this approach, too much water can mix with land areas, and if I narrow the latitude/longitude range on land, I might miss too much land to compute the median meaningfully.
Therefore, with this data, would it be possible to use an IF/ELSE statement to tell R that if the "center point" of each grid cell happens to fall on land, then it would be considered as land (i.e. that would be TRUE - if not, then FALSE)? Even if a grid cell happens to have water mixed with land, but the center point of the grid is on land, that would be considered land. But can R even tell what is land or water in this case?
Thank you, and I would greatly appreciate any assistance!

        [[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
Reply | Threaded
Open this post in threaded view
|

Re: Isolating only land areas on a global map for computing averages

R-sig-geo mailing list
In reply to this post by Tom Philippi
Hi Tom and others,

I was able to use rasterToPolygons(RCP1pctCO2Mean) to convert my raster into polygons, which is now an object called "trans". However, does this command retain the original 138 layers/years of data, but just now in polygon form? 
Also, I am not too clear as to how to generate id values for polygons as row and column indices from this new object. I was looking online for a procedure, but nothing too specific shows how to do this directly. Would it be possible to provide an example in this context? 
Finally, I am unsure of what was meant by "get land polygon SpatialPolygons object "landPoly" with polygons of land only, not ocean." Do you mean the SpatialPolygonsDataframe that was created from rasterToPolygons?
Thank you, and I appreciate your help!
-----Original Message-----
From: Tom Philippi <[hidden email]>
To: rain1290 <[hidden email]>; tephilippi <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Fri, Nov 8, 2019 11:04 pm
Subject: Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages

Rain--Yes, it is possible to do it with your extant raster stack.  In fact, pretty much all reasonable approaches will do that.  Anything you do will create a raster layer with values at every one of your 8192 raster cells: some values will be 0 or NA.
The trivial answer if you only had a small range of latitude, and small raster cell sizes relative to your polygon, would be the raster::rasterize() function to generate a raster mask.
But, with a global raster from +90 to -90, any interpretable answer averaging over area needs to take into account the different area of a grid cell at 70 deg latitude vs 0 deg.  See: https://gis.stackexchange.com/questions/29734/how-to-calculate-area-of-1-x-1-degree-cells-in-a-rasterAt that point, you should go ahead and account for fractions of grid cells over land so your averaging would be over land area.
I'm reticent to give you a complete code solution because without a name, you may well be a student with this as an assignment.  But, my approach would be:
create a SpatialPolygonsDataFrame object "gridPoly" from any of your raster layers via raster::rasterToPolygons()generate an id value for each of those polygons as row & column indices from the raster, or as the cell number.get land polygon SpatialPolygons object "landPoly" with polygons of land only, not ocean.create a new SpatialPolygons object from gIntersect::gridPoly, landPoly, byid = c(TRUE, FALSE))calculate an area of each polygon in that object via geosphere::areaPolygon()  {rgeos::gArea()  only works for projected CRS}create your mask/weights raster layer with either all NA or all 0.either:      parse the id values to row & column values.     use the triplets of row, column, and area to replace the corresponding NA or 0 values in that mask
or:     if you used cell numbers, just use the cell numbers and area values in replacement in that mask
create a second weight raster via raster::area() on one of your raster layers.
raster multiply your polygon-area and your raster::area values to give the actual weighs to use.
This still is an approximation, but likely +/- 1-2%.
If this is still complete gibberish to you, either I need more coffee or you need to consult a good reference on working with spatial data in general.
On Thu, Nov 7, 2019 at 8:25 AM <[hidden email]> wrote:

Hi Tom and others,

Thank you for your response and suggestions! 
Yes, I loaded and used the maptools package previously to create continents on my world map, among other things. I do think that the easiest approach would be to create a raster layer for land, and then water, with the values that I have. However, my precipitation values are globally distributed - every grid cell has a precipitation value for each year (i.e. each grid cell has 138 values/layers/years). So, if I were to create a raster of only land areas, how would I have my grid cells coincide with the land areas only on that raster?
Also, would it be possible to accomplish this with the raster stack that I already have? If so, is there a way to separate all land/water areas this way using the maptools package?
Thanks, again, and I really appreciate your help!
-----Original Message-----
From: Tom Philippi <[hidden email]>
To: rain1290 <[hidden email]>
Cc: r-sig-geo <[hidden email]>
Sent: Thu, Nov 7, 2019 12:44 am
Subject: Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages

The easiest approach would be to create a separate aligned raster layer for land vs water.  There are plenty of coastline polygons available out there (e.g., maptools, rworldmap, rnaturalearth packages): choose one in your raster CRS (or choose one and spTransform() it).  Then, use a grid version of your raster to extract values from that land/water SpatialPolygons object.
1: Your idea of extracting the land/water value at each grid cell centroid gives one estimate.  Instead of TRUE/FALSE, think of the numeric equivalents 1,0,  then using those as weights for averaging across your grid cells.2: A "better" estimate would be to compute the fraction of each grid cell that is land, then use those fractional [0, 1] values as weights to compute a weighted average of precipitation over land.  At 2.8deg grid cells, a lot of heavy rainfall coastal areas would have the grid cell centroid offshore and be omitted by approach #1.3: I recommend that you think hard about averaging across cells in Lat Lon to estimate average precipitation over land.  The actual area of a ~2.8 by 2.8 deg grid cell at the equator is much larger than the same at 70 deg N.  I would spend the extra hour computing the actual area (in km^2) of land in each of your 8192 grid cells, then using those in a raster as weights for whatever calculations you do on the raster stack.  [Or you can cheat, as the area of a grid cell in degrees is a function of only the latitudes, and your required weights are multiplicative.]
Your mileage may vary...
Tom
On Wed, Nov 6, 2019 at 6:18 PM rain1290--- via R-sig-Geo <[hidden email]> wrote:

Hi there,
I am interested in calculating precipitation medians globally. However, I only want to isolate land areas to compute the median. I already first created a raster stack, called "RCP1pctCO2median", which contains the median values. There are 138 layers, with each layer representing one year.  This raster stack has the following attributes:
class       : RasterStack
dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
resolution  : 2.8125, 2.789327  (x, y)
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ...
min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730, 0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687, 0.51857087, 0.41005131, 0.45956529, 0.47497867, ...
max values  :   96.30350,  104.08584,   88.92751,   97.49373,   89.57201,   90.58570,   86.67651,   88.33519,   96.94720,  101.58247,   96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...  

Previously, I was isolating a specific region by specifying a range of longitudes and latitudes to obtain the medians for that region, like this:
expansion1<-expand.grid(103:120, 3:15)lonlataaa<-extract(RCP1pctCO2Median, expansion1)Columnaaa<-colMeans(lonlataaa)

However, with this approach, too much water can mix with land areas, and if I narrow the latitude/longitude range on land, I might miss too much land to compute the median meaningfully.
Therefore, with this data, would it be possible to use an IF/ELSE statement to tell R that if the "center point" of each grid cell happens to fall on land, then it would be considered as land (i.e. that would be TRUE - if not, then FALSE)? Even if a grid cell happens to have water mixed with land, but the center point of the grid is on land, that would be considered land. But can R even tell what is land or water in this case?
Thank you, and I would greatly appreciate any assistance!

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