Shapefile reprojection with "sf" R package doesn't work

Question

I am using the sf R package to work with a shapefile (with a species distribution).

> sp #this is the shapefile

returns:

Simple feature collection with 1 feature and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -11539950 ymin: 4933183 xmax: 11985800 ymax: 8012809
epsg (SRID):    NA
proj4string:    +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
            scntfcN                       geometry
3096 Myodes rutilus MULTIPOLYGON (((5070195 797...

These are the shapefile crs and plot:

> crs(sp) #show the original crs
[1] "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
> plot(sp) #plot the original shapefile

I want to reproject the shapefile from the Mollweide projection to the WGS84 geographic crs:

> sp_reproj <- st_transform(sp, 4326) #reproject the shapefile

> crs(sp_reproj)

returns:

[1] "+proj=longlat +datum=WGS84 +no_defs"

but when I try to plot it, I get this error:

> plot(sp_reproj)
Error in CPL_geos_is_empty(st_geometry(x)) : 
  Evaluation error: IllegalArgumentException: point array must contain 0 or 1 elements.

Furthermore, if I write the shapefile

> st_write(sp_reproj, "Reprojection.shp")

and open it with QGIS, I have this:

Strangely, the holes (lakes) seems to have been filled.

Further notes: 1) the original sp shapefile seems to be valid and not empty

> st_is_valid(sp)
[1] TRUE
> st_is_empty(sp)
[1] FALSE

2) My final aim is not just to plot the reprojected shapefile, but to intersect it with a raster. However, I get the same error as when plotting, so the problem should be in the shapefile not in the operation.

3) Actually, the sp shapefile come from a bigger shapefile with thousand of sf features. This one is the only one that gives me that error.

4) Here the original and reprojected shapefiles

Why I get this error? And why the reprojection fill the holes in the original shapefile?

Answer 1

Your sp object is a MULTIPOLYGON - lets break it into POLYGON objects and see if there's something afoot:

> sp2 = st_cast(sp$geometry,"POLYGON")
> length(sp2)
[1] 40

Its 40 separate POLYGON units. At this point plot(sp2) works and looks fine. Now transform those 40 units and plot:

> sp2t = st_transform(sp2, 4326)
> plot(sp2t)
Error in CPL_geos_is_empty(st_geometry(x)) : 
  Evaluation error: IllegalArgumentException: point array must contain 0 or >1 elements.

Great, we've reproduced your error, which is the first step to fixing it. Now which of the 40 units is the problem? A little testing reveals:

> plot(sp2t[1:39])

works fine. but...

> plot(sp2t[40])
Error in CPL_geos_is_empty(st_geometry(x)) : 
  Evaluation error: IllegalArgumentException: point array must contain 0 or >1 elements.

What is that pesky object?

> sp2t[40]
Geometry set for 1 feature 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 179.9987 ymin: 67.11926 xmax: 179.9987 ymax: 67.11926
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POLYGON ((179.9987 67.11926))

It appears to be a polygon with only one point. What did it start out as?

> sp2[40]
Geometry set for 1 feature 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 9936218 ymin: 7483013 xmax: 10072890 ymax: 7528533
epsg (SRID):    NA
proj4string:    +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
POLYGON ((9936320 7528522, 10072891 7483013, 99...

which is a polygon of four coordinates:

> st_coordinates(sp2[40])
            X       Y L1 L2
[1,]  9936320 7528522  1  1
[2,] 10072891 7483013  1  1
[3,]  9936218 7528533  1  1
[4,]  9936320 7528522  1  1

which transformed becomes one point:

> st_coordinates(sp2t[40])
            X        Y L1 L2
[1,] 179.9987 67.11926  1  1

The original polygon is a tiny sliver at the top right. It has an area that is about one ten-millionth of the total area:

> sum(st_area(sp2))
1.846842e+13 [m^2]
> st_area(sp2[40])
1600110 [m^2]

and it is probably being reprojected to a point because of edge problems with the transformation. Let's try projecting each point in polygon 40:

> p40pts = st_cast(sp2[40],"POINT")
> st_transform(p40pts, 4326)
Geometry set for 4 features  (with 3 geometries empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: 179.9987 ymin: 67.11926 xmax: 179.9987 ymax: 67.11926
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POINT EMPTY
POINT EMPTY
POINT (179.9987 67.11926)
POINT EMPTY

Empty geometries. These are probably right on the edge or even over the edge of that Mollweide projection.

So what to do? I reckon this sliver is negligible, so you can rebuild your spatial data by dropping it and remaking a MULTIPOLYGON object if you need it:

> spfix = sp2t[1:39]
> spfixu = st_union(spfix)
> plot(spfixu, col="brown")

and the lakes are preserved as holes.