3

When projecting to the Azimuthal Equidistant projection in R's rgdal the results seem strange. Take this example:

require(rgdal)
require(maptools)
data(wrld_simpl)
azim_polar = '+proj=aeqd +lat_0=90 +lon_0=0'
azim_orign = '+proj=aeqd +lat_0=0 +lon_0=0'
w_polar = spTransform(wrld_simpl, CRS(azim_polar))
w_orign = spTransform(wrld_simpl, CRS(azim_orign))
par(mai=c(0,0,.2,0), mfrow=c(1, 2))
plot(w_polar, col = 'grey80'); title(azim_polar)
plot(w_orign, col = 'grey80'); title(azim_orign)

enter image description here

Antarctica's polygon should fill outwards to the projection edge, not plot inwards. Similarly for the WGS84 origin plot, shouldn't the outer polygons continue outwards instead of wrapping round? Compare this to the d3.js implementation).

Am I missing something in the process, or is this an rgdal or proj4 issue?

Edit: plotting meridians and parallels to clarify what's going on here

enter image description here

Note that the polar plot is actually correct, but the other is fully wrapped around. I've made lines semi-transparent to show overplotting, and the only lines not overplotting in the second map are the 90° E and W meridians. May be wrong but that suggests to me a problem in the math rather than plotting function..

geotheory
  • 22,624
  • 29
  • 119
  • 196
  • Page 45 of the [Proj.4 documentation](ftp://ftp.remotesensing.org/proj/OF90-284.pdf) is relevant. I've tried those examples, which come out very differently using `rgdal`.. The other relevant resources are [on github](https://github.com/OSGeo/proj.4/wiki). – geotheory Feb 22 '16 at 23:53
  • I would say its a bug in sp's plotting. The coordinates of Antarctica are projected correctly, but sp doesn't know it should really overlay it with the outline circle. What does sp do with polygons that wrap over the horizon with other spherical projections? I suspect it also screws up... – Spacedman Feb 23 '16 at 08:07
  • Hi Barry, yes I think that might explain the polar azimuthal. There's definitely also some erroneous wrapping going on. I've added grid maps above to clarify - note comments. – geotheory Feb 23 '16 at 13:27
  • Looks like this may be due to [a known Proj.4 bug](https://github.com/OSGeo/proj.4/issues/357) that is resolved in `4.9.2`, but `rgdal` (1.0-4, (SVN revision 548)) is still on `4.9.1`.. – geotheory Feb 23 '16 at 15:52
  • The general issue is described here in detail: https://bost.ocks.org/mike/example/ You can fix this sort of stuff in R by re-composing the polygons sensibly, but I don't think there's a general good answer until the topological boundaries are able to be handled, and that's a bigger problem than proj.4 can deal with. – mdsumner Apr 06 '16 at 08:19
  • Maybe able to fix this case by decomposing to small triangles, I'll try that – mdsumner Apr 06 '16 at 08:20

1 Answers1

0

Yes, updating Proj.4 and gdal with homebrew has partly fixed it. Still have a plotting issue with Antarctica however..

enter image description here

geotheory
  • 22,624
  • 29
  • 119
  • 196