5

Using both R and QGIS I got unexpected output from the buffer function. I'm wondering if it has something to do with my coordinate reference system, USA Contiguous Albers Equal Area Conic, because the issue goes away after reprojecting to UTM Zone 15. Can anyone tell me why this would be occurring? Am I missing something basic here?

Distorted buffer output plotted on original Comparing Albers versus UTM

R code for reproducing the error and figures:

library(sp)
library(rgeos)

utm15 = "+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
albers =  "+init=epsg:5070 +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

field_bndry = readWKT("POLYGON ((-95.5176901800000024 42.0119487799999973, -95.5177545500000065 42.0120284899999987, -95.5177719900000000 42.0099468800000011, -95.5177760100000057 42.0097814700000001, -95.5177907599999969 42.0097615399999995, -95.5178162499999956 42.0097575500000033, -95.5178377000000012 42.0097316399999983, -95.5178336799999954 42.0096808199999998, -95.5177719900000000 42.0096439499999974, -95.5177813800000024 42.0094426599999977, -95.5177907599999969 42.0093529800000027, -95.5178149000000047 42.0093350400000034, -95.5179409699999979 42.0092921899999965, -95.5180321600000042 42.0092603000000011, -95.5181086099999987 42.0092383799999993, -95.5181729800000028 42.0092343900000031, -95.5182105300000046 42.0092343900000031, -95.5183714599999973 42.0091546800000017, -95.5186021300000050 42.0090350999999984, -95.5188113499999929 42.0089035600000003, -95.5190527400000065 42.0088158699999994, -95.5192136799999929 42.0086803399999980, -95.5193638800000002 42.0085966400000004, -95.5193638800000002 42.0085488100000006, -95.5178511099999952 42.0091467000000023, -95.5177867400000054 42.0091188000000031, -95.5178725699999944 42.0090829299999982, -95.5178618399999948 42.0090191500000003, -95.5177867400000054 42.0089912500000011, -95.5177760100000057 42.0084372000000030, -95.5186665100000027 42.0084172700000025, -95.5190795700000024 42.0084212499999978, -95.5193102399999958 42.0084332099999997, -95.5194658000000061 42.0084810399999995, -95.5196964699999995 42.0085727199999965, -95.5200880799999936 42.0086404799999968, -95.5202490100000006 42.0086803399999980, -95.5203723900000057 42.0086763600000026, -95.5208283700000038 42.0085049600000033, -95.5221104599999933 42.0084730700000009, -95.5247336599999954 42.0084611099999989, -95.5275392500000038 42.0084571299999965, -95.5275821700000023 42.0089354500000027, -95.5276250800000071 42.0094934800000033, -95.5274963399999990 42.0100674600000019, -95.5272602999999947 42.0103544400000004, -95.5270242699999983 42.0107530299999965, -95.5267882299999940 42.0112472699999984, -95.5265736599999968 42.0115980300000018, -95.5261874199999994 42.0118690600000022, -95.5262517900000034 42.0120125500000015, -95.5226898199999965 42.0120125500000015, -95.5226898199999965 42.0117893400000000, -95.5227327300000013 42.0114704799999998, -95.5229044000000016 42.0112313300000011, -95.5228829400000023 42.0109284099999982, -95.5229687700000056 42.0107051999999968, -95.5231618899999972 42.0104182099999974, -95.5234623000000056 42.0100833999999992, -95.5237841600000053 42.0097964099999999, -95.5239343599999984 42.0094934800000033, -95.5236983300000020 42.0094456500000035, -95.5232477200000005 42.0096210299999981, -95.5230116799999962 42.0100036799999970, -95.5227971099999991 42.0104022700000002, -95.5225396200000034 42.0105298200000021, -95.5223035799999991 42.0105776500000019, -95.5223035799999991 42.0111037799999991, -95.5221533800000060 42.0115342499999969, -95.5221533800000060 42.0120284899999987, -95.5176901800000024 42.0119487799999973))")
proj4string(field_bndry) = CRS("+init=epsg:4326")

field_bndry_nad83 = spTransform(field_bndry, CRS("+init=epsg:4269"))
field_bndry_albers = spTransform(field_bndry_nad83, albers)
field_bndry_utm = spTransform(field_bndry_nad83, utm15)


png("albers_vs_utm.png")
par(mfrow=c(4, 2))
plot(field_bndry_albers, main="Albers", sub="No buffer")
plot(field_bndry_utm, main="UTM 15", sub="No buffer")

plot(
  gBuffer(field_bndry_albers, width = 0),
  sub="0m buffer"
)
plot(
  gBuffer(field_bndry_utm, width=0),
  sub="buffer 0m"
)

plot(
  gBuffer(field_bndry_albers, width = 1),
  sub="1m buffer"
)
plot(
  gBuffer(field_bndry_utm, width=1),
  sub="buffer 1m"
)


plot(
  gBuffer(field_bndry_albers, width = 10),
  sub="10m buffer"
)
plot(
  gBuffer(field_bndry_utm, width=10),
  sub="buffer 10m"
)


dev.off()

png("original_vs_buffer.png")
plot(field_bndry_albers, col="#a6cee3")
plot(gBuffer(field_bndry_albers, width = 30), main="30m buffer", add=T, col="#1f78b4")
legend(
  "topleft",
  c("Original boundary", "30m buffer"),
  fill=c("#a6cee3", "#1f78b4")
)
dev.off()

4

Your polygon isn't valid - it has (at least) a self-intersection in the top right corner:

> gIsValid(field_bndry)
[1] FALSE
Warning message:
In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point -95.517755208090747 42.0119499413616

This means GEOS' idea of what's the inside and what's the outside is unclear, and the buffer doesn't look like it makes sense. You need to fix this geometry before you can proceed.

You can fix it with gMakeValid:

> plot(field_bndry_albers)
> plot(gBuffer(field_bndry_albers,width=30))
> plot(gBuffer(gMakeValid(field_bndry_albers),width=30))

enter image description here

5
  • Wow, I thought I had checked that, my apologies for my lack of due diligence. Thank you! Is there something about the projection to UTM which makes the geometry valid? I'm seeing that gIsValid still throws a warning in the UTM CRS.
    – Dave-Evans
    Apr 25 at 18:32
  • What makes you think the projection to UTM makes the geometry valid? gIsValid(field_bndry_utm) returns FALSE.
    – Spacedman
    Apr 25 at 19:31
  • My apologies, I was unclear. I mean to say what about UTM makes the buffer work as aspected despite the geometry being invalid.
    – Dave-Evans
    Apr 25 at 20:42
  • It looks like maybe it reversed the order of the points in the polygon. Hence it things the little triangle in the top right is now outside the polygon and it seems to disappear in your image instead of being buffered.
    – Spacedman
    Apr 25 at 21:00
  • Wow! Good spot, thank you. I was thinking it must have been something to do with that projection but it seems instead like it's happenstance. Anyhow, lesson learned, when in doubt, use gIsValid. Thanks for your time.
    – Dave-Evans
    Apr 25 at 21:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.