How to draw grid contours within a map

? 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.

--
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