st_transform producing spurious ribbons across region

Question

I am trying to use st_transform on PERZONES shapefile available here . In its native projection, it looks like this:

However, when I run the command:

tansformedPF<-st_transform(PF,"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs") 
##This should be the MODIS sinusoidal projection

I get:

I thought maybe st_transform_proj would fix the problem but it did not. I also tried to plot the transformed file in Earth Engine just in case it was a problem with R's plotting and not the raster itself, but EE also shows the bands. I've used st_transform with the same destination crs with other shapefiles with no problem.

I also tried st_wrap_dateline(PF, options = c("WRAPDATELINE=YES") and st_wrap_dateline(PF, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180") before passing to st_transform, but no luck. The native proj4string is "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs".

Any workarounds to reproject this file without the bands? I am open to another reproject function, but st_transform seems to be the one to use.

Answer 1

Here's the full solution. To split polygons, use st_difference with a very thin buffer - this also makes sure you don't get any edge effects at the tear line.

I'm going to work with a subset of the biggest polygons to save time. Read the data in and choose the big ones:

> perm_full = st_read("./UiO_PEX_PERZONES_5.0_20181128_2000_2016_NH.shp")
> pa = perm_full[st_area(perm_full) > max(st_area(perm_full))/4,]

Next construct a line from the pole up to somewhere near the equator at the top of the projected map. This will be our tear line when we reproject. Make a ten metre wide buffer polygon:

> am = st_linestring(x=rbind(c(0,0),c(0,4000000)))
> amb = st_buffer(am, dist=10)

Now split. When plotted you can see a line which is actually a 10 metre gap sliced through that chunk of what I think is Russia. My apologies to Vladimir Putin.

> pasplit = st_difference(pa$geometry, st_sfc(amb,crs=st_crs(pa)))
> plot(pasplit)

Now reproject using your modis projection...

> modis = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
> plot(st_transform(pasplit, modis))

tested on the full dataset:

Answer 2

Thank you @Spacedman

Alternatively, this also works:

PF %>% 
  st_transform(4326) %>% 
  st_wrap_dateline() %>% 
  st_transform("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs") -> 
  s2

According to edzer here, "The issue arises if the antimeridian is used to split the Earth for projection, as does sinusoidal and also equidistant rectilinear (used by sf for plotting geographic, longlat coordinates). "