how to specify the size of device so that it fits a SpatialPolygonsDataFrame plot exactly

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

how to specify the size of device so that it fits a SpatialPolygonsDataFrame plot exactly

Patrick Giraudoux
Dear listers,

I would like to write jpeg (or png, etc.) maps using the function
graphics:jpeg (png, etc.), but with the device window exactly fitting
the map (= no white strip above, below or on the sides of the map). The
map is derived from a SpatialPolygonsDataFrame (WGS84)

However, even if I set

par(mar=c(0,0,0,0) and pass xaxs='i', yaxs='i' in the plot, there are
still white strips

I tried several tricks but none are fully satisfying.

Trial #1: I computed the ratio of the bounding box and using this ratio
in the function jpeg with the height and width arguments. I did it both
in WGS84 and in UTM47:

ratio<-(max(bbox(ChinaUTM)[1,])-min(bbox(ChinaUTM)[1,]))/(max(bbox(ChinaUTM)[2,])-min(bbox(ChinaUTM)[2,]))
# ratio width/height

then

jpeg("./MapChina.jpeg",height=8,width=8*ratio,res=300,units="in")

par(mar=c(0,0,0,0)) # no margin
plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')

dev.off()

The best is obtained with UTM47 (planar projection), however I get
(almost) rid of any white strip only if I add an additionnal 1.04
coefficient (obtained by trial-error)

Hence:

jpeg("./MapChina.jpeg",height=8,width=8*ratio*1.04,res=300,units="in")

par(mar=c(0,0,0,0)) # no margin
plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')

dev.off()

Trial #2: The other way has been to pick up the parameter values like that:

par(mar=c(0,0,0,0)) # no margin
plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
pt<-par()$pin

ratio<-pt[2]/pt[1]

jpeg("./MapChina.jpeg",height=8,width=8*ratio,res=300,units="in")
par(mar=c(0,0,0,0)) # no margin
plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i') #

dev.off()

Does not work either

Any idea about how to deal with this ? In short how to get the exact
size of a SpatialPolygonDataFrame plot to make the device window exactly
this  size?

Best,

Patrick

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: how to specify the size of device so that it fits a SpatialPolygonsDataFrame plot exactly

Roger Bivand
Administrator
On Sat, 27 Jul 2019, Patrick Giraudoux wrote:

> Dear listers,
>
> I would like to write jpeg (or png, etc.) maps using the function
> graphics:jpeg (png, etc.), but with the device window exactly fitting the map
> (= no white strip above, below or on the sides of the map). The map is
> derived from a SpatialPolygonsDataFrame (WGS84)
>
> However, even if I set
>
> par(mar=c(0,0,0,0) and pass xaxs='i', yaxs='i' in the plot, there are still
> white strips
>
> I tried several tricks but none are fully satisfying.
>
> Trial #1: I computed the ratio of the bounding box and using this ratio in
> the function jpeg with the height and width arguments. I did it both in WGS84
> and in UTM47:
>
> ratio<-(max(bbox(ChinaUTM)[1,])-min(bbox(ChinaUTM)[1,]))/(max(bbox(ChinaUTM)[2,])-min(bbox(ChinaUTM)[2,]))
> # ratio width/height
>
> then
>
> jpeg("./MapChina.jpeg",height=8,width=8*ratio,res=300,units="in")
>
> par(mar=c(0,0,0,0)) # no margin
> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
>
> dev.off()
>
> The best is obtained with UTM47 (planar projection), however I get (almost)
> rid of any white strip only if I add an additionnal 1.04 coefficient
> (obtained by trial-error)
>
> Hence:
>
> jpeg("./MapChina.jpeg",height=8,width=8*ratio*1.04,res=300,units="in")
>
> par(mar=c(0,0,0,0)) # no margin
> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
>
> dev.off()
>
> Trial #2: The other way has been to pick up the parameter values like that:
>
> par(mar=c(0,0,0,0)) # no margin
> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
> pt<-par()$pin
>
> ratio<-pt[2]/pt[1]
>
> jpeg("./MapChina.jpeg",height=8,width=8*ratio,res=300,units="in")
> par(mar=c(0,0,0,0)) # no margin
> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i') #
>
> dev.off()
>
> Does not work either
>
> Any idea about how to deal with this ? In short how to get the exact size of
> a SpatialPolygonDataFrame plot to make the device window exactly this  size?
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?

Best wishes,

Roger

>
> Best,
>
> Patrick
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: 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 :

> On Sat, 27 Jul 2019, Patrick Giraudoux wrote:
>
>> Dear listers,
>>
>> I would like to write jpeg (or png, etc.) maps using the function
>> graphics:jpeg (png, etc.), but with the device window exactly fitting
>> the map (= no white strip above, below or on the sides of the map).
>> The map is derived from a SpatialPolygonsDataFrame (WGS84)
>>
>> However, even if I set
>>
>> par(mar=c(0,0,0,0) and pass xaxs='i', yaxs='i' in the plot, there are
>> still white strips
>>
>> I tried several tricks but none are fully satisfying.
>>
>> Trial #1: I computed the ratio of the bounding box and using this
>> ratio in the function jpeg with the height and width arguments. I did
>> it both in WGS84 and in UTM47:
>>
>> ratio<-(max(bbox(ChinaUTM)[1,])-min(bbox(ChinaUTM)[1,]))/(max(bbox(ChinaUTM)[2,])-min(bbox(ChinaUTM)[2,]))
>> # ratio width/height
>>
>> then
>>
>> jpeg("./MapChina.jpeg",height=8,width=8*ratio,res=300,units="in")
>>
>> par(mar=c(0,0,0,0)) # no margin
>> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
>>
>> dev.off()
>>
>> The best is obtained with UTM47 (planar projection), however I get
>> (almost) rid of any white strip only if I add an additionnal 1.04
>> coefficient (obtained by trial-error)
>>
>> Hence:
>>
>> jpeg("./MapChina.jpeg",height=8,width=8*ratio*1.04,res=300,units="in")
>>
>> par(mar=c(0,0,0,0)) # no margin
>> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
>>
>> dev.off()
>>
>> Trial #2: The other way has been to pick up the parameter values like
>> that:
>>
>> par(mar=c(0,0,0,0)) # no margin
>> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i')
>> pt<-par()$pin
>>
>> ratio<-pt[2]/pt[1]
>>
>> jpeg("./MapChina.jpeg",height=8,width=8*ratio,res=300,units="in")
>> par(mar=c(0,0,0,0)) # no margin
>> plot(China,col="grey",border="white",lwd=2,xaxs='i', yaxs='i') #
>>
>> dev.off()
>>
>> Does not work either
>>
>> Any idea about how to deal with this ? In short how to get the exact
>> size of a SpatialPolygonDataFrame plot to make the device window
>> exactly this  size?
>
> 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?
>
> Best wishes,
>
> Roger


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.

Best,

Patrick

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: how to specify the size of device so that it fits a SpatialPolygonsDataFrame plot exactly

Roger Bivand
Administrator
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: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: 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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo