[R-sig-Geo] Create a simple spatial grid

[Barry Rowlingson]

On Tue, May 20, 2014 at 8:44 PM, Remi Genevest <rgenevest at free.fr> wrote:
> You're right, such a high resolution for all over Europe is unrealistic.
>
> I should rather build a smaller grid for the 20m*20m cells and make an other
> coarser grid for the whole Europe.
>
> But, anyway, how can I create a spatial grid from scratch, by defining CRS,
> resolution and extent of my grid ?
>

 The raster package has functions for doing all that.

# coordinate systems by epsg code - the EU one and lat-long:

eu = CRS("+init=epsg:3035")
ll = CRS("+init=epsg:4326")

# I'm going to make a 101x81 grid of *approximate* 20m squares with
origin at this lat-long point:

origin = SpatialPoints(cbind(-2.783897,54.008399),ll)

# I have to convert that origin point to EU projection coordinates,
which are in metres:

require(rgdal)
origin.eu = spTransform(origin, eu)

# and I'm going to specify my extent. I've rounded the origin
coordinates to the nearest 100:

e = extent(3487700, 3487700+(101*20), 3507300, 3507300 + (81*20))

# now I define a raster based on that extent:

r = raster(e)

# and tell it it has 20m square cells:

res(r)=c(20,20)

# now I can fill it with 81*101 data values:

r[] = runif(81*101)

# set its projection:

projection(r)=eu

# and plot it.

plot(r)

# I can save it as a GeoTIFF:

writeRaster(r,"/tmp/uni.tif")

Now, if I use QGIS to plot that grid over an OpenStreetMap basemap
then it doesn't line up with NS, because the UK is quite far west in
europe and the angular change in a conical projection is extremal at
the west and east limits. If you want NS aligned rasters you need
another CRS, and that will have other distortions in it. C'est la vie.

Barry