Mixing points and sf objects

Dear listers,

Just for curiosity, using one sf object (depts) and one classical
data.frame with two columns bearing coordinates, I wonder why this
sequence works (points and department limits are correctly drawn):

plot(chbrut[,5:6],asp=1)
plot(st_geometry(depts),add=TRUE)

.. but not this one (points of chbrut are not projected - or projected
somewhere else, and are not visible on the map)

plot(st_geometry(depts),reset=FALSE)
points(chbrut[,5:6],pch=19,col=cols)

Best,

Patrick

	[[alternative HTML version deleted]]

On Thu, 3 Mar 2022, Patrick Giraudoux wrote:

Dear listers,

Just for curiosity, using one sf object (depts) and one classical
data.frame with two columns bearing coordinates, I wonder why this
sequence works (points and department limits are correctly drawn):

plot(chbrut[,5:6],asp=1)
plot(st_geometry(depts),add=TRUE)

.. but not this one (points of chbrut are not projected - or projected
somewhere else, and are not visible on the map)

plot(st_geometry(depts),reset=FALSE)
points(chbrut[,5:6],pch=19,col=cols)

Sorry, cannot reproduce:

library(sf)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
nc_xy <- st_transform(nc, "EPSG:32019")
set.seed(1)
pts <- as.data.frame(st_coordinates(st_sample(nc_xy, 100)))

plot(pts, asp=1)
plot(st_geometry(nc_xy), add=TRUE)

and

plot(st_geometry(nc_xy))
points(pts)

seem the same. reset= is irrelevant if the plot method for sfc objects 
is used, I think it applies to sf objects only, when multiple maps may 
be displayed with keys.

Hope this helps,

Roger

Best,

Patrick

????[[alternative HTML version deleted]]

walk  Simple feature collection with 2095 features and 2 fields Geometry 

st_write(walk,dsn="walk.gpx",layer="track_points")  Error: Could not guess driver for walk.gpx

Dear listers,

I am making trials with the sf package trying to write point geometries
into afile. Here is the sf object:

walk  Simple feature collection with 2095 features and 2 fields Geometry

type: POINT Dimension: XY Bounding box: xmin: 5.951268 ymin: 47.29238
xmax: 5.974512 ymax: 47.3015 Geodetic CRS: WGS 84 First 10 features:
track_seg_point_id ele geometry 1 0 243 POINT (5.95348 47.29709) 2 1 243
POINT (5.953438 47.29704) 3 2 244 POINT (5.953219 47.29702) 4 3 243
POINT (5.952884 47.29713) 5 4 243 POINT (5.952492 47.2972) 6 5 243 POINT
(5.952335 47.29743) 7 6 242 POINT (5.951881 47.29732) 8 7 242 POINT
(5.951495 47.29724) 9 8 293 POINT (5.951333 47.29719) 10 9 294 POINT
(5.951388 47.29706)

Writing a shapefile as well as a kml is no problem:

st_write(walk,"walk.shp")

st_write(walk,"walk.kml")

But I cannot go through writing a gpx file, e.g.

st_write(walk,dsn="walk.gpx",layer="track_points")  Error: Could not

guess driver for walk.gpx

I have tried several other combinations all unsuccessful.

Have someone an idea about what goes wrong ?

Best,

Patrick





        [[alternative HTML version deleted]]

walk? Simple feature collection with 2095 features and 2 fields Geometry 

st_write(walk,dsn="walk.gpx",layer="track_points")? Error: Could not guess driver for walk.gpx

walk  Simple feature collection with 2095 features and 3 fields Geometry 

Hi Patrick,

You can add the driver manually.

st_write(walk, dsn = "walk.gpx", layer = "track_points", driver = "GPX")

Hope it helps,
Ahmadou


On Sat, Mar 5, 2022 at 6:53 PM Patrick Giraudoux 
<patrick.giraudoux at univ-fcomte.fr> wrote:

    Dear listers,

    I am making trials with the sf package trying to write point
    geometries
    into afile. Here is the sf object:

    > walk? Simple feature collection with 2095 features and 2 fields

    Geometry
    type: POINT Dimension: XY Bounding box: xmin: 5.951268 ymin: 47.29238
    xmax: 5.974512 ymax: 47.3015 Geodetic CRS: WGS 84 First 10 features:
    track_seg_point_id ele geometry 1 0 243 POINT (5.95348 47.29709) 2
    1 243
    POINT (5.953438 47.29704) 3 2 244 POINT (5.953219 47.29702) 4 3 243
    POINT (5.952884 47.29713) 5 4 243 POINT (5.952492 47.2972) 6 5 243
    POINT
    (5.952335 47.29743) 7 6 242 POINT (5.951881 47.29732) 8 7 242 POINT
    (5.951495 47.29724) 9 8 293 POINT (5.951333 47.29719) 10 9 294 POINT
    (5.951388 47.29706)

    Writing a shapefile as well as a kml is no problem:

    st_write(walk,"walk.shp")

    st_write(walk,"walk.kml")

    But I cannot go through writing a gpx file, e.g.

    > st_write(walk,dsn="walk.gpx",layer="track_points") Error: Could

    not guess driver for walk.gpx

    I have tried several other combinations all unsuccessful.

    Have someone an idea about what goes wrong ?

    Best,

    Patrick





    ? ? ? ? [[alternative HTML version deleted]]

I have made some additionnal trials, with the argument 'driver'
explicit. It appears that the sf object 'walk' must compulsory have two
fields names "track_fid" and "track_seg_id" respectively. Once created :

walk  Simple feature collection with 2095 features and 3 fields Geometry

type: POINT Dimension: XY Bounding box: xmin: 5.951268 ymin: 47.29238
xmax: 5.974512 ymax: 47.3015 Geodetic CRS: WGS 84 First 10 features:
track_fid ele geometry track_seg_id 1 1 243 POINT (5.95348 47.29709) 1 2
1 243 POINT (5.953438 47.29704) 1 3 1 244 POINT (5.953219 47.29702) 1 4
1 243 POINT (5.952884 47.29713) 1 5 1 243 POINT (5.952492 47.2972) 1 6 1
243 POINT (5.952335 47.29743) 1 7 1 242 POINT (5.951881 47.29732) 1 8 1
242 POINT (5.951495 47.29724) 1 9 1 293 POINT (5.951333 47.29719) 1 10 1
294 POINT (5.951388 47.29706) 1


Then,

st_write(walk,dsn="walk.gpx",layer="track_points",driver="GPX")

is written with no error

However, st_write(walk,dsn="walk.gpx",layer="track_points") still gets
'Error: Could not guess driver for walk.gpx'




Le 05/03/2022 ? 20:24, Ahmadou Dicko a ?crit :

Hi Patrick,

You can add the driver manually.

st_write(walk, dsn = "walk.gpx", layer = "track_points", driver = "GPX")

Hope it helps,
Ahmadou


On Sat, Mar 5, 2022 at 6:53 PM Patrick Giraudoux
<patrick.giraudoux at univ-fcomte.fr> wrote:

    Dear listers,

    I am making trials with the sf package trying to write point
    geometries
    into afile. Here is the sf object:

    > walk  Simple feature collection with 2095 features and 2 fields

    Geometry
    type: POINT Dimension: XY Bounding box: xmin: 5.951268 ymin: 47.29238
    xmax: 5.974512 ymax: 47.3015 Geodetic CRS: WGS 84 First 10 features:
    track_seg_point_id ele geometry 1 0 243 POINT (5.95348 47.29709) 2
    1 243
    POINT (5.953438 47.29704) 3 2 244 POINT (5.953219 47.29702) 4 3 243
    POINT (5.952884 47.29713) 5 4 243 POINT (5.952492 47.2972) 6 5 243
    POINT
    (5.952335 47.29743) 7 6 242 POINT (5.951881 47.29732) 8 7 242 POINT
    (5.951495 47.29724) 9 8 293 POINT (5.951333 47.29719) 10 9 294 POINT
    (5.951388 47.29706)

    Writing a shapefile as well as a kml is no problem:

    st_write(walk,"walk.shp")

    st_write(walk,"walk.kml")

    But I cannot go through writing a gpx file, e.g.

    > st_write(walk,dsn="walk.gpx",layer="track_points") Error: Could

    not guess driver for walk.gpx

    I have tried several other combinations all unsuccessful.

    Have someone an idea about what goes wrong ?

    Best,

    Patrick





            [[alternative HTML version deleted]]

        [[alternative HTML version deleted]]

Hi,

I have a trouble with the combination of layout and plot.stars. e.g.

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

plot(mypoly,col=mypal(12),border=mypal(12))

I expect that the first two plots display in region #1, the second added 
to the first, and the third plot in region #2. However, this is not what 
happens: actually, the third plot displays in region 1 erasing the 
others. I understand that plot.stars when not "added" does not respect 
the layout definition (and displays its own regions), and that my 
problem comes from the way plot.stars deals with that.

Has anyone an idea about a workaround ?

On Tue, 15 Mar 2022, Patrick Giraudoux wrote:

 Hi,

 I have a trouble with the combination of layout and plot.stars. e.g.

 nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
 layout.show(nf)

 plot(st_geometry(mydept),col="grey",border="grey90")
 plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

You have already noticed that sf and stars, like raster and terra, modify the 
assumptions of base plot methods, as graphics::filled.contour(), unless some 
ordering and argument conditions are met, crucially the non-base reset= 
argument. I do not think that you can use layout() at all.

library(stars)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
bir74_rast <- st_rasterize(nc["BIR74"])
plot(bir74_rast, reset=FALSE)
plot(st_geometry(nc),border="grey90", add=TRUE)
gridGraphics::grid.echo()
library(grid)
g <- grid.grab()
plot(nc["BIR74"])
gridGraphics::grid.echo()
gv <- grid.grab()
grid.newpage()
gridExtra::grid.arrange(g, gv, ncol=2)

Grabbing the base graphics device state lets you use 
gridExtra::grid.arrange() to place multiple graphics objects; here I haven't 
tried to constrain aspect or relative sizes. I don't think that the plot 
methods in sf and stars play well with layout, because they use it themselves 
internally.

Hope this helps,

Roger

 plot(mypoly,col=mypal(12),border=mypal(12))

 I expect that the first two plots display in region #1, the second added
 to the first, and the third plot in region #2. However, this is not what
 happens: actually, the third plot displays in region 1 erasing the others.
 I understand that plot.stars when not "added" does not respect the layout
 definition (and displays its own regions), and that my problem comes from
 the way plot.stars deals with that.

 Has anyone an idea about a workaround ?

On Tue, 15 Mar 2022, Patrick Giraudoux wrote:

Hi,

I have a trouble with the combination of layout and plot.stars. e.g.

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

You have already noticed that sf and stars, like raster and terra, 
modify the assumptions of base plot methods, as 
graphics::filled.contour(), unless some ordering and argument 
conditions are met, crucially the non-base reset= argument. I do not 
think that you can use layout() at all.

library(stars)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
bir74_rast <- st_rasterize(nc["BIR74"])
plot(bir74_rast, reset=FALSE)
plot(st_geometry(nc),border="grey90", add=TRUE)
gridGraphics::grid.echo()
library(grid)
g <- grid.grab()
plot(nc["BIR74"])
gridGraphics::grid.echo()
gv <- grid.grab()
grid.newpage()
gridExtra::grid.arrange(g, gv, ncol=2)

Grabbing the base graphics device state lets you use 
gridExtra::grid.arrange() to place multiple graphics objects; here I 
haven't tried to constrain aspect or relative sizes. I don't think 
that the plot methods in sf and stars play well with layout, because 
they use it themselves internally.

Hope this helps,

Roger

plot(mypoly,col=mypal(12),border=mypal(12))

I expect that the first two plots display in region #1, the second 
added to the first, and the third plot in region #2. However, this is 
not what happens: actually, the third plot displays in region 1 
erasing the others. I understand that plot.stars when not "added" 
does not respect the layout definition (and displays its own 
regions), and that my problem comes from the way plot.stars deals 
with that.

Has anyone an idea about a workaround ?

Great ! Thanks Roger.? On this basis, I have a way to explore the issue
now. Will give a feed-back on the list once done.
Best,
Patrick

Le 15/03/2022 ? 11:27, Roger Bivand a ?crit?:

On Tue, 15 Mar 2022, Patrick Giraudoux wrote:

Hi,

I have a trouble with the combination of layout and plot.stars. e.g.

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

You have already noticed that sf and stars, like raster and terra,
modify the assumptions of base plot methods, as
graphics::filled.contour(), unless some ordering and argument
conditions are met, crucially the non-base reset= argument. I do not
think that you can use layout() at all.

library(stars)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
bir74_rast <- st_rasterize(nc["BIR74"])
plot(bir74_rast, reset=FALSE)
plot(st_geometry(nc),border="grey90", add=TRUE)
gridGraphics::grid.echo()
library(grid)
g <- grid.grab()
plot(nc["BIR74"])
gridGraphics::grid.echo()
gv <- grid.grab()
grid.newpage()
gridExtra::grid.arrange(g, gv, ncol=2)

Grabbing the base graphics device state lets you use
gridExtra::grid.arrange() to place multiple graphics objects; here I
haven't tried to constrain aspect or relative sizes. I don't think
that the plot methods in sf and stars play well with layout, because
they use it themselves internally.

Hope this helps,

Roger

plot(mypoly,col=mypal(12),border=mypal(12))

I expect that the first two plots display in region #1, the second
added to the first, and the third plot in region #2. However, this is
not what happens: actually, the third plot displays in region 1
erasing the others. I understand that plot.stars when not "added"
does not respect the layout definition (and displays its own
regions), and that my problem comes from the way plot.stars deals
with that.

Has anyone an idea about a workaround ?

	[[alternative HTML version deleted]]

In addition to what Roger wrote, you can use e.g.

layout(matrix(1:4,2))
library(stars)
s = st_as_stars(L7_ETMs)
image(s[,,,1])
plot(s[,,,1], key.pos = NULL, reset = FALSE)
plot(s[,,,1], key.pos = NULL, reset = FALSE, main = NULL)

to fill the sub-plots of layout() incrementally. Note the key.pos and 
reset arguments to plot.stars(): they make sure plot.stars doesn't 
mess with the layout settings.

On 15/03/2022 11:45, Patrick Giraudoux wrote:

Great ! Thanks Roger.? On this basis, I have a way to explore the issue
now. Will give a feed-back on the list once done.
Best,
Patrick

Le 15/03/2022 ? 11:27, Roger Bivand a ?crit?:

On Tue, 15 Mar 2022, Patrick Giraudoux wrote:

Hi,

I have a trouble with the combination of layout and plot.stars. e.g.

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

You have already noticed that sf and stars, like raster and terra,
modify the assumptions of base plot methods, as
graphics::filled.contour(), unless some ordering and argument
conditions are met, crucially the non-base reset= argument. I do not
think that you can use layout() at all.

library(stars)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
bir74_rast <- st_rasterize(nc["BIR74"])
plot(bir74_rast, reset=FALSE)
plot(st_geometry(nc),border="grey90", add=TRUE)
gridGraphics::grid.echo()
library(grid)
g <- grid.grab()
plot(nc["BIR74"])
gridGraphics::grid.echo()
gv <- grid.grab()
grid.newpage()
gridExtra::grid.arrange(g, gv, ncol=2)

Grabbing the base graphics device state lets you use
gridExtra::grid.arrange() to place multiple graphics objects; here I
haven't tried to constrain aspect or relative sizes. I don't think
that the plot methods in sf and stars play well with layout, because
they use it themselves internally.

Hope this helps,

Roger

plot(mypoly,col=mypal(12),border=mypal(12))

I expect that the first two plots display in region #1, the second
added to the first, and the third plot in region #2. However, this is
not what happens: actually, the third plot displays in region 1
erasing the others. I understand that plot.stars when not "added"
does not respect the layout definition (and displays its own
regions), and that my problem comes from the way plot.stars deals
with that.

Has anyone an idea about a workaround ?

????[[alternative HTML version deleted]]

In addition to what Roger wrote, you can use e.g.

layout(matrix(1:4,2))
library(stars)
s = st_as_stars(L7_ETMs)
image(s[,,,1])
plot(s[,,,1], key.pos = NULL, reset = FALSE)
plot(s[,,,1], key.pos = NULL, reset = FALSE, main = NULL)

to fill the sub-plots of layout() incrementally. Note the key.pos and 
reset arguments to plot.stars(): they make sure plot.stars doesn't 
mess with the layout settings.

On 15/03/2022 11:45, Patrick Giraudoux wrote:

Great ! Thanks Roger.? On this basis, I have a way to explore the issue
now. Will give a feed-back on the list once done.
Best,
Patrick

Le 15/03/2022 ? 11:27, Roger Bivand a ?crit?:

On Tue, 15 Mar 2022, Patrick Giraudoux wrote:

Hi,

I have a trouble with the combination of layout and plot.stars. e.g.

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

You have already noticed that sf and stars, like raster and terra,
modify the assumptions of base plot methods, as
graphics::filled.contour(), unless some ordering and argument
conditions are met, crucially the non-base reset= argument. I do not
think that you can use layout() at all.

library(stars)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
bir74_rast <- st_rasterize(nc["BIR74"])
plot(bir74_rast, reset=FALSE)
plot(st_geometry(nc),border="grey90", add=TRUE)
gridGraphics::grid.echo()
library(grid)
g <- grid.grab()
plot(nc["BIR74"])
gridGraphics::grid.echo()
gv <- grid.grab()
grid.newpage()
gridExtra::grid.arrange(g, gv, ncol=2)

Grabbing the base graphics device state lets you use
gridExtra::grid.arrange() to place multiple graphics objects; here I
haven't tried to constrain aspect or relative sizes. I don't think
that the plot methods in sf and stars play well with layout, because
they use it themselves internally.

Hope this helps,

Roger

plot(mypoly,col=mypal(12),border=mypal(12))

I expect that the first two plots display in region #1, the second
added to the first, and the third plot in region #2. However, this is
not what happens: actually, the third plot displays in region 1
erasing the others. I understand that plot.stars when not "added"
does not respect the layout definition (and displays its own
regions), and that my problem comes from the way plot.stars deals
with that.

Has anyone an idea about a workaround ?

????[[alternative HTML version deleted]]

Coming back to the example which initiated the issue, I have found the 
trick following your hints Edzer, but it makes a strange rational. Here 
we go:

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_geometry(fzone1),col=mypal(12)[1],border=NA,add=TRUE)
plot(st_geometry(fzone2),col=mypal(12)[2],border=NA,add=TRUE)
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",key.pos=NULL,reset=FALSE,add=TRUE)
plot(st_geometry(mydept),col=NA,border="grey90",add=TRUE)

mybox<-bbox2sf(n=0.5,s=0,w=0,e=10,crs=2154) # crs is just crap here, I 
do not need it thenafter
mypoly<-st_make_grid(mybox,n=c(12,1),what="polygons")
plot(mypoly,col=mypal(12),border=mypal(12))
plot(mybox,add=TRUE)

Makes exactly what I was intending to make. The critical point was to 
explicitely keep reset=FALSE in the 4th plot (actually a plot.stars). If 
not, the 6th plot is plotted in region #1 erasing the previous plots...

Some personal remarks:
- in the first plot, I ommited the argument reset=FALSE, however this 
does not make a problem (maybe because I just plot a geometry ?)

- the 4th plot is definitely strange with reset=FALSE and add=TRUE 
together, isn't it ?

Thanks Roger and Ezder for bailing me out once again... :-)
Patrick



Le 15/03/2022 ? 16:02, Edzer Pebesma a ?crit?:

In addition to what Roger wrote, you can use e.g.

layout(matrix(1:4,2))
library(stars)
s = st_as_stars(L7_ETMs)
image(s[,,,1])
plot(s[,,,1], key.pos = NULL, reset = FALSE)
plot(s[,,,1], key.pos = NULL, reset = FALSE, main = NULL)

to fill the sub-plots of layout() incrementally. Note the key.pos and 
reset arguments to plot.stars(): they make sure plot.stars doesn't 
mess with the layout settings.

On 15/03/2022 11:45, Patrick Giraudoux wrote:

Great ! Thanks Roger.? On this basis, I have a way to explore the issue
now. Will give a feed-back on the list once done.
Best,
Patrick

Le 15/03/2022 ? 11:27, Roger Bivand a ?crit?:

On Tue, 15 Mar 2022, Patrick Giraudoux wrote:

Hi,

I have a trouble with the combination of layout and plot.stars. e.g.

nf <- layout(matrix(c(1,2),2,1,byrow = TRUE), c(3,3), c(3,1), TRUE)
layout.show(nf)

plot(st_geometry(mydept),col="grey",border="grey90")
plot(st_rasterize(ztot),col=mypal(12),main="",breaks="equal",add=TRUE)

You have already noticed that sf and stars, like raster and terra,
modify the assumptions of base plot methods, as
graphics::filled.contour(), unless some ordering and argument
conditions are met, crucially the non-base reset= argument. I do not
think that you can use layout() at all.

library(stars)
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
bir74_rast <- st_rasterize(nc["BIR74"])
plot(bir74_rast, reset=FALSE)
plot(st_geometry(nc),border="grey90", add=TRUE)
gridGraphics::grid.echo()
library(grid)
g <- grid.grab()
plot(nc["BIR74"])
gridGraphics::grid.echo()
gv <- grid.grab()
grid.newpage()
gridExtra::grid.arrange(g, gv, ncol=2)

Grabbing the base graphics device state lets you use
gridExtra::grid.arrange() to place multiple graphics objects; here I
haven't tried to constrain aspect or relative sizes. I don't think
that the plot methods in sf and stars play well with layout, because
they use it themselves internally.

Hope this helps,

Roger

plot(mypoly,col=mypal(12),border=mypal(12))

I expect that the first two plots display in region #1, the second
added to the first, and the third plot in region #2. However, this is
not what happens: actually, the third plot displays in region 1
erasing the others. I understand that plot.stars when not "added"
does not respect the layout definition (and displays its own
regions), and that my problem comes from the way plot.stars deals
with that.

Has anyone an idea about a workaround ?

????[[alternative HTML version deleted]]