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
how to specify the size of device so that it fits a SpatialPolygonsDataFrame plot exactly
6 messages · Roger Bivand, Patrick Giraudoux
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 R-sig-Geo at r-project.org 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: Roger.Bivand at nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
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
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 ?
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()
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