Error in creating map objects

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: https://rpubs.com/gergness/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:

https://international.ipums.org/international/gis_yrspecific_1st.shtml

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)
)

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.

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 https://www.jla-data.net/eng/lines-from-points/ 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.

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!

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

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.

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!

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

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.