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

[Patrick Giraudoux]

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




	[[alternative HTML version deleted]]