Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages
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
To: rain1290 ; tephilippi
Cc: r-sig-geo
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 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
To: rain1290
Cc: r-sig-geo
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 aver
Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages
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
To: rain1290 ; tephilippi
Cc: r-sig-geo
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 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
To: rain1290
Cc: r-sig-geo
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
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-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 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
> To: rain1290
> Cc: r-sig-geo
> 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 ch
Re: [R-sig-Geo] Isolating only land areas on a global map for computing averages
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 To: rain1290 Cc: r-sig-geo 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 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 mai
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@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@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
