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

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: enter image description here

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.

  • Sorry, the link to the shapefile is here: doi.pangaea.de/10.1594/PANGAEA.888600?format=html#download – webbe May 14 at 12:53
  • that is just some of your polygons going "behind" the map, clip to the valid region of your new projection – Ian Turton May 14 at 12:58
  • @Ian Turton if I understand you correctly, I should just clip to the land regions ("valid region"?). However, that doesn't solve the problem of the bands across Russia. – webbe May 14 at 13:03
  • Any polygon that crosses a "tear" (discontinuity) in a reprojection needs to be split into multiple polygons precisely along the tear line. (sorry no time for full answer at this point) – Spacedman May 14 at 13:50

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)

enter image description here

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))

enter image description here

tested on the full dataset:

enter image description here


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") -> 

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). "

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.