Skip to content

Create a simple spatial grid

2 messages · Barry Rowlingson, Remi Genevest

#
On Tue, May 20, 2014 at 8:44 PM, Remi Genevest <rgenevest at free.fr> wrote:
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
#
Thanks a lot Barry. Your post perfectly matches what I was looking for.

Well now if I want to do a coarser grid covering Europe, what I should do is
doing more or less the same thing, just to change the resolution, right? 

res(r)=c(100,100)  # an approximate 100meters squared cells

Then I got my raster in a EU coordinates system, and I would like to
transform it into a WGS 84 (so as adding other layers would become easier).
I tried this way, but it doesn't work :

# change the CRS into WGS84 (ll)
r.new<-projectRaster(r,crs="+init=epsg:4326") 
# adapt the extent to Europe
ext = extent(-10.417,31.917,34.083,71.083)
extent(r.new)<-ext
# plot the map
plot(r.new,axes=TRUE)
# Here I get a plot that is not projected in a proper way.

So, I wondered if avoiding the projectRaster() function would be alright :

# Change only the extent adapt the extent to Europe
ext = extent(-10.417,31.917,34.083,71.083)
extent(r)<-ext
projection(r)<-"+init=epsg:4326"
plot(r,axes=TRUE)

Am I wrong ?!

Best regards
R?mi



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Create-a-simple-spatial-grid-tp7586481p7586490.html
Sent from the R-sig-geo mailing list archive at Nabble.com.