[R-sig-Geo] How to draw grid contours within a map
Lyndon Estes
lestes at princeton.edu
Tue Feb 15 15:43:26 CET 2011
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
More information about the R-sig-Geo
mailing list