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 <tephili...@gmail.com> To: rain1290 <rain1...@aim.com>; tephilippi <tephili...@gmail.com> Cc: r-sig-geo <r-sig-geo@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 <rain1...@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 <tephili...@gmail.com> To: rain1290 <rain1...@aim.com> Cc: r-sig-geo <r-sig-geo@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@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
