[R-sig-Geo] How to draw grid contours within a map
Lyndon Estes
lestes at princeton.edu
Tue Feb 15 17:48:09 CET 2011
Hi Thiago,
It seems to work just fine with geographic coordinates as well. I am
working at the moment with a GCS WGS84 shape of Africa (af.gsc), and
the following seems to work fine:
# Africa shapefile download from
# http://www.internationalmapping.com/index/Goodies_Africa.html
p.base <- "/path/to/Africa.shp
setwd(p.base)
gcs.proj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
af <- readOGR(dsn = "African_Coastline.shp", layer = "African_Coastline")
af.gsc <- spTransform(af, gcs.proj) # Transform it GSC
ex <- bbox(af.gsc)
# 10 degree grid
gr <- createPolyGrid(x = ex[1, 1], y = ex[2, 1], h = ex[2, 2] - ex[2,
1], w = ex[1, 2] - ex[1, 1], size = 10)
plot(af.gsc)
plot(gr, add = T)
You might want to check the accuracy of the grids (caveat emptor), or
perhaps a wiser person than myself can comment on the soundness of
this approach.
Cheers, Lyndon
On Tue, Feb 15, 2011 at 11:31 AM, Thiago Veloso <thi_veloso at yahoo.com.br> wrote:
> Hi Lyndon,
>
> Thank you very much for sharing your function and for the tips as well. The code is very useful, but I'll have to make slight modifications in order to properly plot my data...
>
> The grids I need to draw span the following coordinates:
>
>>fcst$lon[121:125]
>>[1] -60.0 -57.5 -55.0 -52.5 -50.0
>>fcstn$lat[23:27]
>>[1] -35.0 -32.5 -30.0 -27.5 -25.0
>
> One solution could be adapting your function to use degrees instead of metres. By doing this, the arguments would be something like:
>
> fcstgrid <- createPolyGrid(x = -60, y = -35, h = 10, w = 10, size = 2.5)
>
> "GridTopology" and related help pages don't mention how to use degrees instead of metres (or my R immaturity prevented me from finding)...
>
> Any tips to help me working on that adaptation?
>
> Best regards,
>
> Thiago.
>
>
> --- On Tue, 15/2/11, Lyndon Estes <lestes at princeton.edu> wrote:
>
>> From: Lyndon Estes <lestes at princeton.edu>
>> Subject: Re: [R-sig-Geo] How to draw grid contours within a map
>> To: "Thiago Veloso" <thi_veloso at yahoo.com.br>
>> Cc: "R SIG list" <r-sig-geo at stat.math.ethz.ch>
>> Date: Tuesday, 15 February, 2011, 12:43
>> Hi Thiago,
>>
>> I was doing something a bit similar to what you are doing
>> just
>> recently. Instead of creating a raster of the grids,
>> I created
>> polygon grids, which is then quite simple to overlay on
>> other
>> polygons.
>>
>> Here's a function I used to create the polygon grids:
>>
>> createPolyGrid <- function(x, y, h, w, size, prj) {
>> # Creates a polygon grid in meters
>> # Args:
>> # x: X coordinate of lower left corner of
>> lower left of polygon, in meters
>> # y: Y coordinate of lower left corner of
>> lower left of polygon, in meters
>> # h: The desired North-South extent of the
>> grid
>> # w: The desired East-West extent of the
>> grid
>> # size: Grid polygon cell size in meters
>> # prj: Object holding the crs string for
>> the desired projection (optional)
>> # Returns:
>> # A SpatialPolygonsDataFrame with square
>> grids, with North-South and
>> East-West distances rounded to
>> # avoid fractionating grid cells. An
>> attribute data frame of 1
>> column holding polygon identifiers is created
>>
>> if(missing(x)) {
>> stop("missing x")
>> } else if(missing(y)) {
>> stop("missing y")
>> } else if(missing(h)) {
>> stop("missing h")
>> } else if(missing(w)) {
>> stop("missing w")
>> } else if(missing(size)) {
>> stop("missing size")
>> }
>>
>> # Define numbers of row and columns in grid,
>> rounding to nearest whole numbers
>> xrow <- round(w / size, 0)
>> yrow <- round(h / size, 0)
>>
>> # Define coordinates for center of lower left cell
>> xcen <- x + size / 2
>> ycen <- y + size / 2
>>
>> # Create grid topology
>> grd <- GridTopology(c(xcen, ycen), c(size, size),
>> c(xrow, yrow))
>>
>> # Turn it into polygons
>> if(missing(prj)) {
>> grd.shp <-
>> as.SpatialPolygons.GridTopology(grd)
>> } else if(!missing(prj)) {
>> grd.shp <-
>> as.SpatialPolygons.GridTopology(grd, prj)
>> }
>>
>> # Create data frame and then SpatialPolygons into
>> SpatialPolygonsDataFrame
>> grd.dat <-
>> as.data.frame(matrix(1:length(grd.shp), ncol = 1))
>> colnames(grd.dat) <- "polID"
>> row.names(grd.dat) <-
>> as.character(grd.dat$polID)
>> grd.out.shp <- SpatialPolygonsDataFrame(grd.shp,
>> grd.dat, match.ID = FALSE)
>>
>> return(grd.out.shp)
>> }
>>
>>
>> # Here's an example of how to make the grid with the
>> function above
>> and plot it over another polygon, based on
>> # the meuse dataset (meuse.sr example taken from
>> http://r-spatial.sourceforge.net/gallery/)
>>
>> library(gstat)
>>
>> data(meuse.riv)
>> meuse.sr =
>> SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse.riv")))
>> # Polygon of river
>>
>> # Find extent of river shape and use the coordinates to
>> create polygon grid
>> ex <- bbox(meuse.sr)
>> gr <- createPolyGrid(x = ex[1, 1], y = ex[2, 1], h =
>> ex[2, 2] - ex[2,
>> 1], w = ex[1, 2] - ex[1, 1], size = 500)
>>
>> # Plot
>> plot(meuse.sr)
>> plot(gr, add = T)
>>
>> Hope this helps.
>>
>> Cheers, Lyndon
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> On Tue, Feb 15, 2011 at 7:22 AM, Thiago Veloso <thi_veloso at yahoo.com.br>
>> wrote:
>> > Dear R-Siggers,
>> >
>> > In order to illustrate the spatial domain of the
>> General Circulation Model (GCM) output data I am using, I
>> need to display some grids within a map. They need to be
>> 2.5°x2.5° lat-lon over specific coordinates. Actually, I
>> am trying to partially reproduce Figure 1 inside Challinor
>> et al., 2005 (image attached).
>> >
>> > To load and display the shapefile I use the
>> following code:
>> >
>> > library(maptools)
>> >
>> southern=readShapePoly("D:/Maps/papers-in-progress/regiao_sul.shp")
>> > plot(southern)
>> >
>> > Following, I create the "empty" grids to draw:
>> >
>> > library(raster)
>> > rprecip<-raster(xmn=-60,xmx=-50,ymn=-35,ymx=-25)
>> > res(rprecip)<-2.5
>> >
>> > But I am stuck on the task of displaying those
>> grids on the map.
>> >
>> > Could anybody please point any directions?
>> >
>> > Best,
>> >
>> > Thiago.
>> >
>> >
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo at r-project.org
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>> >
>>
>>
>>
>> --
>> Lyndon Estes
>> Research Associate
>> Woodrow Wilson School
>> Princeton University
>> +1-609-258-2392 (o)
>> +1-609-258-6082 (f)
>> +1-202-431-0496 (m)
>> lestes at princeton.edu
>>
>
>
>
>
>
>
>
--
Lyndon Estes
Research Associate
Woodrow Wilson School
Princeton University
+1-609-258-2392 (o)
+1-609-258-6082 (f)
+1-202-431-0496 (m)
lestes at princeton.edu
More information about the R-sig-Geo
mailing list