Banding appears when reprojecting shapefiles using R?

Question

I have downloaded a shape file (Marine Ecoregions (last one on the page)) from the nature conservancy http://maps.tnc.org/gis_data.html

I am then using R and the sf package to project this shapefile. If I plot this map on NAD83, WGS84 etc., the map looks good. If I project it to Eckert4, Lambert Azithmuthal, or Mollweide I get banding across the map.

Any idea why this happens?

Here is some reproducible code:

You will need to download the shapefile from the link above and add them to same folder as the function below

# set chooseProjection to one of the following" "WGS84","Eckert4","molli","lambertaz"
project_maps <- function(chooseProjection) {

  projection <- switch(chooseProjection,
                       "WGS84" = "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0",
                       "NAD83" = "+init=epsg:4269 +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0",
                       "NAD27" = "+init=epsg:4267 +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat",
                       "Eckert4" = "+proj=eck4",
                       "molli" = "+proj=moll",
                       "tranmerc"= "+proj=tmerc",
                       "lambertaz"="+proj=laea +x_0=0 +y_0=0 +lon_0=0 +lat_0=0")
  MEco <- sf::st_read(dsn=here::here(),layer="meow_ecos",quiet=T)
  MEco <- sf::st_transform(MEco,crs=projection)

  # plot ecoregion 
  plotMEco <- ggplot2::ggplot() + 
    ggplot2::geom_sf(data = MEco, ggplot2::aes(fill = ECOREGION)) + 
    ggplot2::theme(legend.position = "none") +
    ggplot2::ggtitle("Marine Ecoregions")
  print(plotMEco)
}

Sample figures:

If i crop the bounding box using

MEco <- sf::st_crop(MEco,c(ymin = -89.99, ymax = 89.99, xmin=-179.99, xmax=179.99))

Then the Eckert4 projection displays just fine. However the Lambert Azithmuthal is still strange.

possible duplicate of https://gis.stackexchange.com/questions/169667/fix-distortions-when-reprojecting-cshapes-world-maps — Ian Turton, Jul 02 '19 at 14:58
@IanTurton I have looked at that thread. While the issue is similar, it is not at the poles but at the equator in one projection (Eckert4) and not sure with the other (Lambert azithmuthal). Maps included — andy, Jul 02 '19 at 16:06
looks like the same issue - just different data, you have marine zones that cross the anti-meridian. — Ian Turton, Jul 02 '19 at 16:17
@IanTurton. The output from sf object has the bounding box set at xmin ymin xmax ymax -180.0000 -89.9000 180.0000 86.9194. This already seems to be within the limits. what am i missing? — andy, Jul 02 '19 at 16:22
Its not the antimeridian, its anything that gets stretched madly by a projection. In your Mollweide example, its feature number 143. The long chord across the top left of the globe is the result of a small line segment of polygon 143 being reprojected. — Spacedman, Jul 02 '19 at 16:25
@Spacedman A few questions: how did you troubleshoot to identify 143 as the issue? How can i rectify the problem? Is this a result of a poorly created shapefile or is the shapefile just fine and that this result is just what happens under certain projections? I don't fully understand (clearly ... ;)) — andy, Jul 02 '19 at 17:11
Some projections are unstable at their edges. I suspect that the features have coordinates exactly at 180/-180. Either the projection algorithm is sometimes flipping sides (-180 goes to +20 million or whatever) or the coordinates are slightly smaller than -180 or larger than +180 like 180.000000000012. — mkennedy, Jul 02 '19 at 17:18
I found 143 by a binary search. Something like plot(ecor$geometry[a:b]) then narrow a and b down by halves. Might look at possible fixes tomorrow. — Spacedman, Jul 02 '19 at 18:17

score 2 · Accepted Answer · answered Jul 03 '19 at 14:10

2

The problem with this:

plot(st_transform(MEco$geom,laea))

is that some of the vertices are quite far apart in lat-long, and even more far apart in LAEA, so you get a big straight line joining them rather than following the true projected shape. Add more vertices to your shapes with st_segmentize:

plot(st_transform(st_segmentize(MEco$geom,10000),laea))

and the problem disappears as the more finely vertexed polygon transforms correctly. The parameter value (10000 above) needs setting according to your data and projection. Too large and you don't get any effect, too small and you add zillions of vertices that don't help.

answered Jul 03 '19 at 14:10

Spacedman

63,755
5
81
115

This is great, thanks. So what are best practices in creating shapefiles? Should the creator worry about these potential projection problems or is this kind of thing the norm and the user is expected to make corrections to suit their projection? cheers – andy Jul 03 '19 at 15:43
Not sure the creator can be expected to deal with all the possible projections - look at some of the non-continuous "interrupted" ones, or the Dymaxion. Ideally the software doing the projection should split polygons or recognise your situation but that's not done at present... – Spacedman Jul 03 '19 at 16:27

Banding appears when reprojecting shapefiles using R?

1 Answers1