4
votes

i would like to create a buffer around a line-shapefile with wgs84 coordinates.

I prepared a shapefile with a single line segment and Datum: D_WGS_1984. After that, i loaded the .shp into R using the 'readOGR' command.

After that i tried the gBuffer method out of the rgeos-package for computing the buffer:

gBuffer(l2, width=1.0, quadsegs=5, capStyle="ROUND", joinStyle="ROUND", mitreLimit=0.01)) 
Warning:
In gBuffer(l2, width = 1, quadsegs = 5, capStyle = "ROUND", joinStyle = "ROUND",  :
Spatial object is not projected; GEOS expects planar coordinates

Obviously the command has a problem with the coordinates. I tried some approaches, but didn't find a solution.

Another example i found for a buffer around points was the following, but i'm not sure how to use it in my case: http://r-sig-geo.2731867.n2.nabble.com/compute-buffer-from-point-shapefile-to-have-shapefile-td4574666.html

Any ideas?

Best regards, Stefan

//update:

Reduced to the relevant part, here is the code:

require("rgeos")
require("rgdal")

l2=readOGR(dsn="C:/Maps", layer="osm_ms")

proj4string(l2) <- CRS("+proj=longlat")
l2.trans <- spTransform(l2, CRS("+proj=longlat"))
summary(l2.trans)

> Object of class SpatialLinesDataFrame
> Coordinates:
>   min       max
> x  7.478942  7.772171
> y 51.840318 52.058856
> Is projected: FALSE 
> proj4string : [+proj=longlat +ellps=WGS84]
> Data attributes:

plot(l2.trans)
plot(gBuffer(l2.trans, width=1.0, quadsegs=5, capStyle="ROUND", joinStyle="ROUND", mitreLimit=0.01))

Probably the line:

Is projected: FALSE is the reason for the problem, but i'm not sure how to use the spTranform and how to find the correct projection.

2
Show us what you tried, right we do not know what you already did, other than that you did something.Paul Hiemstra

2 Answers

7
votes

Think about your buffer size of 1.0 units. That would be, in lat-long, 1.0 degrees. Which doesn't make much sense, since 1 degree N-S is different from 1 degree E-W. GEOS is trying to stop you from doing something it thinks is a bit odd.

So you can fool it by assigning it almost any projected coordinate CRS string, doing the buffer, then assigning it back. The underlying numbers wont change. For example if you do:

 proj4string(l2) = CRS("+init=epsg:27700")

you are lying to the system that the numbers are Uk Grid metres. Then you do the buffer, and you give the units which you know are really degrees. GEOS does the calculation using the numbers, assuming the points are on a grid (and not a sphere). Then you just set the CRS back. The numbers don't change.

Actually, there appears to be a proper ESPG code for projected lat-long, so ignore the previous line, and do:

 proj4string(l2) = CRS("+init=epsg:3395") 

EPSG code database here: http://www.epsg-registry.org/ - note the "scope" in the details for epsg:3395 is 'very small scale mapping'.

If you really want a buffer in metres, then you have to convert to a metric projection (for work in the UK I always use epsg:27700) using spTransform which does change the underlying numbers.

1
votes

The answer to your question is in the error you get. Your data needs to be projected, not in a latlon system. Take a look at spTransform from the rgdal package to perform the projection.