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