Error in creating map objects

[lambort]

Hi everyone!

So I'm new to R, but I'm trying to map in/out provincial migration with R Studio. I've been using the code provided here: RPubs - Using ipumsr. The code works well, but when I get to the part where I am creating the map objects, I get this message when I run prov_point:

Warning message:
In st_point_on_surface.sf(.) :
st_point_on_surface assumes attributes are constant over geometries of x

And so I subsequently get this message when I try to run prov_lines:

Error in Ops.sfc(xi, xj) : operation > not supported

At this point, I've tried looking at the other forums to see why this is going wrong, but I can't figure it out. Anyone have suggestions?

Here's the code:

Create map objects ----

China 2000 Shape File from:

IPUMS International

vignette("ipums-geography", package = "ipumsr")

proj4_china <- "+proj=lcc +lat_1=18 +lat_2=24 +lat_0=21 +lon_0=114 +x_0=500000 +y_0=500000 +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs "
proj4_latlon <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

provinces <- read_ipums_sf("geo1_cn2000.zip", verbose = FALSE) %>%
ms_simplify(keep = 0.01)

#This is where I get the warning message
prov_point <- provinces %>%
st_transform(proj4_china) %>%
st_point_on_surface() %>%
st_transform(proj4_latlon) %>%
select(IPUM2000, ADMIN_NAME)

#This is where I get the message that the operation is not supported
prov_lines <- prov_point %>%
st_transform(proj4_latlon) %>%
crossing(., .) %>%
as_data_frame() %>%
rename(GEO1START = IPUM2000, NAMESTART = ADMIN_NAME, GEO1END = IPUM20001, NAMEEND = ADMIN_NAME1) %>%
filter(GEO1START != GEO1END) %>%
mutate(
geometry = gcIntermediate(as(geometry, "Spatial"), as(geometry1, "Spatial")) %>%
map(st_linestring) %>%
st_sfc()
) %>%
select(-geometry1) %>%
st_as_sf()

here is where I can no longer run the code without prov_lines
prov_lines_data <- ipums_shape_inner_join(
migration_by_province,
prov_lines,
by = c("GEO1START", "GEO1END"),
verbose = FALSE
) %>%
mutate(
ANON = N_IN < 20 | N_OUT < 20,
WT_N_IN = ifelse(ANON, NA, WT_N_IN),
WT_N_OUT = ifelse(ANON, NA, WT_N_OUT),
LINE_WIDTH = sqrt(abs(WT_N_NET)) / 100,
COLOR = ifelse(WT_N_NET < 0, "#f1a340", "#998ec3"),
LABEL = map(paste0(
NAMESTART, " - ", NAMEEND,
"
In: ", comma(WT_N_IN),
"
Out: ", comma(WT_N_OUT),
"
Net: ", comma(WT_N_NET)
), HTML)
)

[jlacko]

It is not an issue with you going wrong, but of old code - the sf::st_point_on_surface() warning is a non issue; the real problem occurs in creation of the prov_lines object - the author of your script uses tidyr::crossing() on two {sf} data frames (or rather uses one data frame to cross with itself to generate all combinations). Current version of {tidyr} no longer works for this use case (the geometry column gets in the way).

This use case - lines connecting centers of all regions, for future highlighting - is a common one, and you should be able to find a working implementation somewhere.
Unfortunately I am not able to point you to one right now though.

[jlacko]

I did not have much success in locating a nice walkthrough for the points to lines use case. As the problem is an interesting, and common, one I ended up writing one myself.

I this blog post Lines from Points · Jindra Lacko I propose a function that can be used as an one liner in your script to generate the data frame of province ids & labels.

As the Chinese provinces use case is a rather specific one I am demonstrating the use case on the well known shapefile of North Carolina counties; it ships with the {sf} package, and as a consequence is widely available. The technique should be easily transferable though.

[lambort]

Thank you so much! This is extremely helpful, I'm going to try out this out and see how it goes. I had tried searching for another example and just couldn't find a good one, so I'm so thankful that you ended up writing one!

[jlacko]

Thank you; I was in a bit of a rush when I wrote the my first reply, but when I returned to the subject later I was also surprised by the scarcity of good examples for this (as I thought...) simple and common use case

The closest is https://docs.ropensci.org/stplanr/ but it was not quite what I had in mind...

So go ahead, give it a try - and if you fail give me a shout. I am thinking of writing a small package around the topic, so user feedback is most welcome.

[lambort]

My map is a bit of a mess (working on that), but it worked! Even I could figure out how to use the function you typed up in about a half hour (still super new to this whole thing). Thank you again!

[jlacko]

This is not a mess - this is exactly what I would expect given the data. You are doing great!
A little finetuning will be required, but you are on the right track. I'm glad to hear that