[R-sig-Geo] Isolating only land areas on a global map for computing averages

[r@i@1290 m@iii@g oii @im@com]

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 <tephilippi using gmail.com>
To: rain1290 <rain1290 using aim.com>; tephilippi <tephilippi using gmail.com>
Cc: r-sig-geo <r-sig-geo using r-project.org>
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 <rain1290 using aim.com> 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 <tephilippi using gmail.com>
To: rain1290 <rain1290 using aim.com>
Cc: r-sig-geo <r-sig-geo using r-project.org>
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 <r-sig-geo using r-project.org> 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
R-sig-Geo using r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



	[[alternative HTML version deleted]]