Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to compute a Line-Buffer with SpatialLinesDataFrame

Tags:

r

gis

spatial

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.

like image 592
nomeus Avatar asked Feb 22 '23 03:02

nomeus


1 Answers

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.

like image 96
Spacedman Avatar answered Feb 24 '23 12:02

Spacedman