Error in plotting with sf_centroid

sbl_bah · February 2, 2020, 8:38pm

Hello,

You will find a warning that I have when I use the sf package.
This warning 2 imply a big problem : when I want to plot (X,Y) points in my background.
these points are not plotted on the shape centroid (they have the same crs)
See snapshot

Have you an idea?

> library(dplyr)
> library(sf)
> 
> map<-st_read("map.gpkg")
Reading layer `countries' from data source `.../map.gpkg' using driver `GPKG'
Simple feature collection with 232 features and 16 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -55.52139 xmax: 180 ymax: 83.59961
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
> 
> pt_centroid<-map %>% sf::st_centroid()
Warning messages:
1: In st_centroid.sf(.) :
  st_centroid assumes attributes are constant over geometries of x
2: In st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon) :
  st_centroid does not give correct centroids for longitude/latitude data
> pts<-sf::st_coordinates(pt_centroid)

jlacko · February 2, 2020, 9:04pm

Without a reproducible example this will be a sort of a guesswork; however I suspect there is an issue with projection of your data. See the hint st_centroid() is giving you in the form of a warning.

The map.gpkg seems to be in WGS84 (EPSG:4326) and that is a geometric CRS; for a centroid to work you need projected data. You need to sf::st_transform() to a different projection - which one will depend on your data, with which we are not familiar with, but pick one with eastings and northings.

To double check consider running st_is_longlat() on your object - as long as it returns TRUE the centroid will not work.

williaml · February 2, 2020, 9:17pm

Is that meant to be Australia on its side?

Perhaps transform it to a projected geographic system.
Something like st_read("map.gpkg") %>% st_transform(3577) which is GDA94 / Australian Albers.

sbl_bah · February 3, 2020, 4:14pm

I take another background where i have the same problem.
I use a function flowmap that has a graphical output. The graphic output (lines) is created from the background map with the sf functions. The lines of the graph connect the centroid
When I want to overlay the background with the graphic the centroids are not aligned.
I can't get the same reference frame to get a good overlay. I think the problem is a graphical setting I don't know.

You can see that the nodes of the chart should be centered on the centroids of the background.

par(bg="NA")
plot(st_geometry(map),col="blue")   # background
plot(df.point$geometry,col="yellow",pich=16,add=TRUE) # centroid of the background
par(new=TRUE)
flowmap(tab=flows,  # Function that's created the line with the centroid of the map  
        fij="Fij",origin.f = "i",destination.f = "j",
        bkg= map,code="EPT_NUM",nodes.X="X",nodes.Y = "Y",
        filter=FALSE)

dev.off()

sbl_bah · February 3, 2020, 4:34pm

Thank you for your answer.

The offset (apparently in X) does not seem to be due to a problem of crs, nor to the use of a local projection (reprojection did not change anything).

It seems rather that it is related to the margins, because when I change the parameter by(mar=c(0,0,1,0)), the offset appears more important.

Maybe it is necessary to add a parameter on the margins?
I can't figure out which one.
Do you have an idea?

williaml · February 3, 2020, 8:09pm

Is this a publicly available dataset? It would be good to have a look at it.

sbl_bah · February 5, 2020, 6:30am

Hello,
It's possible to send you background and dataset.

So, i would like to know if you have understood the description of my problem.

sbl_bah · February 8, 2020, 1:12pm

Hello,

You will find hereafter the reprex() of my function test, shape file and picture of the result.
I would like that the point green and red must be in the same position to ensure the overlay of the segments and background map.
Note : I tried to upload a shape file but it's not possible.
Thanks in advance to your help!

rm(list=ls())
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3

#data(flowdata)
map <- st_read("MGP_TER.shp")
#> Cannot open data source MGP_TER.shp
#> Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : Open failed.

nodes <- data.frame(stringsAsFactors = FALSE,
                    Xi= c(651786.1,651786.1,651786.1),
                    Yi=c(6862044,6862044,6862044),
                    Xj=c(651786.1,662465.4,665350.7),
                    Yj=c(6862044,6858149,6851330)
)
p.X <- sf::st_as_sf(x = nodes, coords = c("Xi", "Yi"), crs = 4326)
p.Y<- sf::st_as_sf(x = nodes, coords = c("Xj", "Yj"), crs = 4326)

nodes<-cbind(p.Y,p.X)  

test<-function(tab){
  plot(tab$geometry, col = "green", lwd = 2)
  plot(tab$geometry.1, add=TRUE ,col = "green", lwd = 2)
  trace <- arrows(tab$Xi, tab$Yi, tab$Xj, tab$Yj, code = 0, col = "black")
}

dev.off()
#> null device 
#>           1

par(bg = "NA")
plot(st_geometry(map), col = "blue") # background
#> Error in st_geometry(map): objet 'map' introuvable
plot(st_centroid(map), add=TRUE,col = "red")
#> Error in st_centroid(map): objet 'map' introuvable
par(new=TRUE)
#> Warning in par(new = TRUE): appel de par(new=TRUE) sans graphe
test(nodes)
#> Warning in st_is_longlat(x): bounding box has potentially an invalid value
#> range for longlat data

sbl_bah · February 9, 2020, 7:05am

Hi, have you an Idea about my problem?
Thanks in advance

williaml · February 9, 2020, 10:53am

Is MGP_TER.shp publicly available?

Did you notice the errors when you used st_read? Maybe that is the issue.

sbl_bah · February 9, 2020, 12:26pm

You can find hereafter the link of MGP_TER.shp here

Any error when I used st_read!

williaml · February 9, 2020, 9:42pm

I assume that you are just trying to reproduce this from the cartograflow package. You should have mentioned that in your post.

Have you looked at the vignette?
https://cran.csiro.au/web/packages/cartograflow/vignettes/cartograflow.html
or the one on the GitHub page? GitHub - fbahoken/cartogRaflow: Tips and tricks for flow and movement mapping

sbl_bah · February 10, 2020, 6:18am

Hello,
I know-well cartograflow i'm author.
I want to upgrade cartograflow to pass to sf instead of sp.
And also i want to modify thé output of thé flowmap with only the graph without background map, it's my issue!

williaml · February 10, 2020, 8:17am

It looks like a good package.

williaml · February 10, 2020, 8:36am

Are you sure it is not something wrong with the crs in the points?

p.X <- sf::st_as_sf(x = nodes, coords = c("Xi", "Yi"), crs = 4326)

You have it set to EPSG:4326 but the they should be in a projected coordinate system. You are looking to convert them to 4326 but you haven't defined what crs they were initially (they are not in 4326).

> p.X 
Simple feature collection with 3 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 651786.1 ymin: 6862044 xmax: 651786.1 ymax: 6862044
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Whereas for 4326, they should be something like this (different dataset):

Simple feature collection with 15 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 144.3703 ymin: -38.30431 xmax: 146.5315 ymax: -37.70646
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

jlacko · February 10, 2020, 10:00am

I concur with your opinion; the bbox of this object is not consistent with EPSG:4326 - you would expect latitude to range between -90 and +90, and longitude between -180, +180.

> p.X 
Simple feature collection with 3 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 651786.1 ymin: 6862044 xmax: 651786.1 ymax: 6862044
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

sbl_bah · February 12, 2020, 7:10pm

thank you for answers.
I have just done the exercise again, projecting the two entities (areas and points) in the same projection system.
The problem does not come from a no coincidence of the CRS. The superposition works perfectly well outside my "test" function:

plot(st_geometry(map),col="blue")
plot(st_geometry(nodes),col="green",add=TRUE)
arrows(nodes$Xi, nodes$Yi, nodes$Xj, nodes$Yj, code = 1, col = "black")

output is great:

When I use my test function, to generate lines from the same nodes, the overlay no longer works. It's as if the areas and the nodes are no longer in the same projection system... Which is not the case.

test<-function(tab){
  plot(st_geometry(tab), col = "green", lwd = 2)
  #plot(tab$geometry.1, add=TRUE ,col = "green", lwd = 2)
  arrows(tab$Xi, tab$Yi, tab$Xj, tab$Yj, code = 0, col = "black")
}
dev.off()
par(bg = "NA")
plot(st_geometry(map),col="red")
plot(st_centroid(map),col="black",add=TRUE)
par(new=TRUE)
seg<-test(tab=nodes)

output is false:

It seems that the overlay of entities in sf doesn't work for entities that come out of a function, ie I can't overlay the same nodes on the zones anymore.

What do you think about this?

jlacko · February 12, 2020, 8:15pm

Now this is weird. The final shape looks plausible, but seems sorta kinda rescaled compared to the original one.

I have seen an unexpected behaviour like this on a workshop recently - it was the result of adding to a plot a manually created sfg geometry, just for the heck of it - and did not pay much attention to it then.

I don't have the time to follow this through right now, but I am intrigued... I will have a look at this later.

sbl_bah · February 15, 2020, 8:18am

Hi,
Is it possible for you to get some références of the workshop? I think it will be interesting to read about this subject.
Thanks in advance

sbl_bah · February 17, 2020, 6:37pm

Hi,
I found the solution to my problem.

In test function i added 'plot(sf :: st_geometry()...)' and After that i can overlay the map background.

Thanks for your help