1

I'd like to have a map of land above 50N in the stereographic projection. I use the following code but as you can see, some land is cut out.

  proj<- "+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 +ellps=WGS84
    +towgs84=0,0,0" 


    data("wrld_simpl", package = "maptools")                                                                            
    wm <- crop(wrld_simpl, extent(-180, 180, 50, 90)) 
    wm <- spTransform(wm, CRSobj = CRS(proj))  

    ### plot land
    plot(wm, col="grey", border="grey")

#add wrapper circle
pts=SpatialPoints(rbind(c(-180,50),c(0,50),c(180,80),c(180,80)), crs(data))
my_line = SpatialLines(list(Lines(Line(cbind(seq(-180,180,0.5),rep(50,721))), ID="outer")), CRS("+init=epsg:4326"))
my_line1<-st_as_sf(my_line)
bounds<-st_transform(my_line1,crs =proj)
plot(bounds,col="gray42",add=TRUE)

enter image description here

1 Answer 1

2

What's happening is that the crop line is a straight line in lat-long with few vertices. spTransform only transforms vertices, so a straight line between two points in lat-long becomes a straight-line in projected coordinates, which is what you see.

To fix, you need to densify straight lines with few vertices into lines with more vertices, and you can use the smoothr package function densify which does a pretty good job here. Install with install.packages("smoothr") as usual and library(smoothr) to get it.

> wm <- crop(wrld_simpl, extent(-180, 180, 50, 90)) 
> plot(spTransform(densify(wm),proj))

enter image description here

1
  • Works like a charm. Thank you, @Spacedman!
    – webbe
    Nov 25, 2019 at 1:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.