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

[Patrick Giraudoux]

Le 29/07/2019 à 11:12, Roger Bivand a écrit :
> 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 


Yes, indeed. Thanks. When one understands fully what's happens, the 
easier to adapt...

Now I can go ahead cleanly...

Best,

Patrick