[R-sig-Geo] how to specify the size of device so that it fits a SpatialPolygonsDataFrame plot exactly

[Roger Bivand]

On Sun, 28 Jul 2019, Patrick Giraudoux wrote:

>>  Le 28/07/2019 à 17:08, Roger Bivand a écrit : This feels like a section in
>>  ASDAR, ch. 3, and code chunks 24-27 in
>>  https://asdar-book.org/book2ed/vis_mod.R. Could you please try that first
>>  by way of something reproducible?
>
> Le 28/07/2019 à 17:15, Patrick Giraudoux a écrit :
>>
>>  Ups... How  can I have missed this chapter of my bible ;-) ? (must admit I
>>  had been on google first). Will re-read it carefully and come back to the
>>  list with a solution or a reproducible example, indeed.
>
>
> OK. Read it again, I was not totally lost. Here is a reproducible example. To 
> ease reproducibility with simple objects, I will use two bounding boxes.  
> bbChina in WGS84, bbChinaUTM47 in UTM47. I want a window fitting the WGS84, 
> and you'll see I get it through a strange way.
>
> bbChina <- new("SpatialPolygons", polygons = list(new("Polygons", Polygons = 
> list( new("Polygon", labpt = c(104.8, 35.95), area = 2372.28, hole = FALSE, 
> ringDir = 1L, coords = structure(c(73, 73, 136.6, 136.6, 73, 17.3, 54.6, 
> 54.6, 17.3, 17.3), .Dim = c(5L, 2L)))), plotOrder = 1L, labpt = c(104.8, 
> 35.95), ID = "1", area = 2372.28)), plotOrder = 1L, bbox = structure(c(73, 
> 17.3, 136.6, 54.6), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"), c("min", 
> "max"))), proj4string = new("CRS", projargs = "+init=epsg:4326 +proj=longlat 
> +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
>
> bbChinaUTM47 <- new("SpatialPolygons", polygons = list(new("Polygons", 
> Polygons = list( new("Polygon", labpt = c(856106.391943348, 4317264.60126758 
> ), area = 30651262771540.1, hole = FALSE, ringDir = 1L, coords = 
> structure(c(-2331000.09677063, -2331000.09677063, 4043212.88065733, 
> 4043212.88065733, -2331000.09677063, 1912947.1678777, 6721582.03465746, 
> 6721582.03465746, 1912947.1678777, 1912947.1678777), .Dim = c(5L, 2L)))), 
> plotOrder = 1L, labpt = c(856106.391943348, 4317264.60126758), ID = "1", area 
> = 30651262771540.1)), plotOrder = 1L, bbox = structure(c(-2331000.09677063, 
> 1912947.1678777, 4043212.88065733, 6721582.03465746), .Dim = c(2L, 2L), 
> .Dimnames = list(c("x", "y"), c("min", "max"))), proj4string = new("CRS", 
> projargs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 
> +towgs84=0,0,0"))
>
> Then let's go:
>
> Example #1: here, being straightforward we get two indesirable white strips 
> on the sides:
>
> width1<-max(bbox(bbChina)[1,])-min(bbox(bbChina)[1,])
> height1<-max(bbox(bbChina)[2,])-min(bbox(bbChina)[2,])
> ratio<-width1/height1
> ratio
> windows(height=8,width=8*ratio)
> par(mar=c(0,0,0,0))
> plot(bbChina,col="grey",xaxs='i', yaxs='i')
> dev.off()
>
> Example #2: computing the ratio on UTM47, but plotting WGS84 (strange), I get 
> a better fit but with still two small white strips up and down.
>
> width1<-max(bbox(bbChinaUTM47)[1,])-min(bbox(bbChinaUTM47)[1,])
> height1<-max(bbox(bbChinaUTM47)[2,])-min(bbox(bbChinaUTM47)[2,])
> ratio<-width1/height1
> ratio
> windows(height=8,width=8*ratio)
> par(mar=c(0,0,0,0))
> plot(bbChina,col="grey",xaxs='i', yaxs='i') # no data range extention (xaxs 
> and yaxs parameter)
>
> dev.off()
>
> Example #3: multiplying the ratio by 1.04, I get a good fit
>
> width1<-max(bbox(bbChinaUTM47)[1,])-min(bbox(bbChinaUTM47)[1,])
> height1<-max(bbox(bbChinaUTM47)[2,])-min(bbox(bbChinaUTM47)[2,])
> ratio<-width1/height1
> ratio
> windows(height=8,width=8*ratio*1.04)
> par(mar=c(0,0,0,0))
> plot(bbChina,col="grey",xaxs='i', yaxs='i')
>
> dev.off()
>
> Looks like the issue has something to do with the way CRS are handled when 
> plotting objects, mmmh ? Tricky isn't it ?
>

Yes, and the section in the book only discusses projected objects, as 
geographical coordinates are stretched N-S proportionally to the distance 
from the Equator. For the UTM47 object, I have:

library(sp)
bbChinaUTM47 <- 
SpatialPolygons(list(Polygons(list(Polygon(matrix(c(-2331000.09677063, 
-2331000.09677063, 4043212.88065733, 4043212.88065733, -2331000.09677063, 
1912947.1678777, 6721582.03465746, 6721582.03465746, 1912947.1678777, 
1912947.1678777), ncol=2))), ID="1")), proj4string=CRS("+proj=utm 
+zone=47")) # you had longlat, so triggering the stretching.

x11() # or equivalent
dxy <- apply(bbox(bbChinaUTM47), 1, diff)
dxy
ratio <- dxy[1]/dxy[2]
ratio
pin <- par("pin")
pin
par(pin=c(ratio * pin[2], pin[2]), xaxs="i", yaxs="i")
plot(bbChinaUTM47)
box()

where the box overlaps the SP object. To finesse:

c(ratio * pin[2], pin[2])
dev.off()
X11(width=6.85, height=5.2)
par(mar=c(0,0,0,0)+0.1)
pin <- par("pin")
par(pin=c(ratio * pin[2], pin[2]), xaxs="i", yaxs="i")
plot(bbChinaUTM47)
box()
dev.off()

From plot.Spatial(), asp is set to 1/cos((mean(ylim) * pi)/180 for 
geographical coordinates, where ylim is a possibly modified version of the 
N-S bounding box. This makes it harder to automate, as you'd need to 
manipulate dxy[2] above to match. So for projected objects, the book 
approach works, for non-projected objects you'd need an extra step.

Hope this helps,

Roger

>
>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en