[R-sig-Geo] How to draw grid contours within a map

[Lyndon Estes]

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