NA handling in `st_interpolate_aw`

[jsta]

If some of the attributes of the "from" polygons in st_interpolate_aw are NA it seems there is no way to avoid them being propagated to the st_interpolate_aw output. Does it make sense to add an na.rm = TRUE argument?

nc = st_read(system.file("shape/nc.shp", package="sf"))
g = st_make_grid(nc, n = c(20,10))
a2 = st_interpolate_aw(nc["BIR74"], g, extensive = TRUE)

nc[sample(seq_len(nrow(nc)), 5), "BIR74"] <- NA
a3 = st_interpolate_aw(nc["BIR74"], g, extensive = TRUE)

plot(a2["BIR74"])
plot(nc["BIR74"])
plot(a3["BIR74"])

[edzer]

Would (should) that be equivalent to removing the corresponding feature altogether?

[jsta]

Hmm, yes it has the same effect. Maybe this issue will be a signpost for future users.

[mountainMath]

I would agree that na.rm=TRUE would be useful. Take the situation where several attributes are to be aggregated, only some of which are NA in some instances. Right now I pipe the first argument through replace(., is.na(.), 0) every time I call st_interpolate_aw, but I feel like a na.rm=TRUE options would be cleaner and more transparent.

[leahbevis]

Hi @mountainMath -- could you tell me how you "pipe the first argument through..." to address this situation? Can you provide sample code? I think am having exactly the same problem: I am aggregating values from numerous smaller polygons that overlap larger polygons, using st_interpolate_aw. Frequently the data that I'm aggregating is missing for 1 of the multiple small polygons, and in those cases I'd like to take the mean of the non-missing values and return that mean to the larger polygon. You have a solution for this? I'm sure what you said above is enough for a seasoned R / dplyr user but I don't use dplyr / piping much.

[leahbevis]

Hey @mountainMath, sorry a follow-up. I realized that -- while I still don't know how to execute the syntax -- your solution is to replace NA with 0, which I guess works when you are taking a sum (extensive = TRUE) but not when you are taking a mean (extensive = FALSE). I don't suppose you or anyone else has a solution for ignoring NA values when taking a mean? I.e., the mean over only the non-NA overlapping polygon values?

[edzer]

Could you pls check & report back whether with this fix, passing na.rm=TRUE works?

[leahbevis]

I re-installed the sf package, and ran this code:

### Demo on st_interpolate_aw
nc = st_read(system.file("shape/nc.shp", package="sf"))
g = st_make_grid(nc, n = c(20,10))

ggplot(data = nc) +
  geom_sf(aes(fill = AREA)) +
  geom_sf(data = g, lwd = .5, color = "#000000",fill = NA)
  
nc$AREA[c(30, 45)] <- NA
nc$AREA[c(30, 46)] <- NA
nc$AREA[c(30, 47)] <- NA
nc$AREA[c(30, 48)] <- NA

ggplot(data = nc) +
  geom_sf(aes(fill = AREA)) +
  geom_sf(data = g, lwd = .5, color = "#000000",fill = NA)

a1 = st_interpolate_aw(nc["AREA"], g, extensive = FALSE)
a2 = st_interpolate_aw(nc["AREA"], g, extensive = FALSE, na.rm=TRUE)

ggplot(data = a1) +
  geom_sf(aes(fill = AREA)) +
  geom_sf(data = g, lwd = .5, color = "#000000",fill = NA)

ggplot(data = a2) +
  geom_sf(aes(fill = AREA)) +
  geom_sf(data = g, lwd = .5, color = "#000000",fill = NA)

The interpolation to g still returns NA for any g polygon that overlaps 1 or more NA polygons in NC. (I.e., those 2 plots above are identical.) So if all I should have to do is re-install sf using install.packages("sf"), then I think I didn't work. But perhaps I'm missing something. I'm working on R version 4.2.3 (2023-03-15) if relevant.

[edzer]

Pls try the dev version, using

install.packages('sf', repos = c('https://r-spatial.r-universe.dev', 'https://cloud.r-project.org'))

[leahbevis]

Ok, I did that just now, used that code above and restarted R and Rstudio for good measure, and that a1 and a2 objects are still identical. I did get this warning when I installed the package: Warning in install.packages : installation of package ‘sf’ had non-zero exit status. Does this mean the install didn't go through properly?

[leahbevis]

Here's an example of the problem with not having na.rm=TRUE pass through. In this pic here you can see the data (gray borders) that I am interpolating up to the black borders. In this pic here you can see the interpolation. Part of the problem may be that the 1970 (large black border) shapefile is less high quality, and the borders aren't matching up properly with the 2020 (gray border) shapefile. Nothing to be done about that, as far as I know. So even though there is little missing data to begin with, the interpolated object ends up with a LOT of missing data.

[edzer]

With this change, I get this for a2:

[edzer]

But t.b.h., I'm still not sure this is the right thing to do.

[leahbevis]

Ugh. I have now updated R, and that sample code works. So, great, in a way. However, due to your phase out of RGDAL, I now can no longe read in shapefiles with readOGR, and I'm not entirely sure how long it's going to take to transition the code that I need for a major presentation on Friday. I was doing a whole set of manipulation of a shapefile using sp, and I imagine it would take me ~10 hours to recreate this code in sf, hours I absolutely do not have. So I can't even GET to to the point of analysis now, as I won't even be able to load and process my data.

Is there any way to re-load the old version of R, 4.2.3? For a mac? Maybe from here?

I know this is NOT your problem, and I really do appreciate the original help. But coincidentally this is just the worst possible timing for this to happen. If you can advise me on un-doing this update I would be so so so incredibly grateful.