Add statistics functions

[jzemmels]

Still to-do:

• unpack "values" and "percentiles" column when computation_type includes "percentile"
• function argument descriptions -- do we want to generalize the get_params, get_description, and base_url functions to pull the docs for /statistics too?
• arg type checking (is that done in other waterdata functions?)
• pull apart get_statistics_data into separate functions to match other read_waterdata patterns
• Figure out if other read_ogc_data helpers are needed:
  • rejigger_cols
  • deal_with_empty
  • switch_arg_id
  • switch_properties_id
• resolve merge conflicts
• (Nice to have) determine if there's a more efficient way to join everything together than the final return_list <- do.call(rbind, return_list_tmp) as it's a slight bottleneck.

And there's Laura's list of to-dos here: https://code.usgs.gov/water/dataRetrieval/-/issues/450

[jzemmels]

To the discussion of whether base R vs. data.table is faster: I've implemented a data.table version of the get_statistics_data function and temporarily put in the read_waterdata_stats_datatable.R script until we decide if we want to add data.table as a dependency in this PR.

I've profiled the rbind-y base R version of the function against the data.table version and the comparison is noteworthy. The data.table version is much faster and uses way-way less memory.

Here's the benchmark code I ran:

bench::mark(
  "data.table" = {
    dat1 <- 
      get_statistics_data_data_table(args = list(approval_status = "approved", 
                                            # state_code = "US:42",
                                            county_code = "US:42:103",
                                            parameter_code = "00095", 
                                            mime_type = "application/json",
                                            computation_type = c("arithmetic_mean", "percentile")),
                                service = "Intervals")
    
    list(
      non_missing_ave = mean(dat1$value, na.rm = TRUE),
      missing_obs = sum(is.na(mean(dat1$value)))
    )
  },
  "base R" = {
    dat2 <-
  get_statistics_data(args = list(approval_status = "approved", 
                                  # state_code = "US:42",
                                  county_code = "US:42:103",
                                  parameter_code = "00095", 
                                  mime_type = "application/json",
                                  computation_type = c("arithmetic_mean", "percentile")),
                      service = "Intervals")
    
    list(
      non_missing_ave = mean(dat2$value, na.rm = TRUE),
      missing_obs = sum(is.na(dat2$value))
    )
  }, min_iterations = 5
)

[ldecicco-USGS]

To the discussion of whether base R vs. data.table is faster: I've implemented a data.table version of the get_statistics_data function and temporarily put in the read_waterdata_stats_datatable.R script until we decide if we want to add data.table as a dependency in this PR.

I've profiled the rbind-y base R version of the function against the data.table version and the comparison is noteworthy. The data.table version is much faster and uses way-way less memory.

Here's the benchmark code I ran:

bench::mark(
  "data.table" = {
    dat1 <- 
      get_statistics_data_data_table(args = list(approval_status = "approved", 
                                            # state_code = "US:42",
                                            county_code = "US:42:103",
                                            parameter_code = "00095", 
                                            mime_type = "application/json",
                                            computation_type = c("arithmetic_mean", "percentile")),
                                service = "Intervals")
    
    list(
      non_missing_ave = mean(dat1$value, na.rm = TRUE),
      missing_obs = sum(is.na(mean(dat1$value)))
    )
  },
  "base R" = {
    dat2 <-
  get_statistics_data(args = list(approval_status = "approved", 
                                  # state_code = "US:42",
                                  county_code = "US:42:103",
                                  parameter_code = "00095", 
                                  mime_type = "application/json",
                                  computation_type = c("arithmetic_mean", "percentile")),
                      service = "Intervals")
    
    list(
      non_missing_ave = mean(dat2$value, na.rm = TRUE),
      missing_obs = sum(is.na(dat2$value))
    )
  }, min_iterations = 5
)

No surprises there, if it wasn't for RDB I would have switched many years ago. Let's go for it. We can start swapping out other code later (off the top of my head I can think of a several places where I'd been tempted in the past). Let's leave those types of edits for a dedicated PR later though.

[ldecicco-USGS]

Let's make sure combined gets returned as a data frame. There are some subtle differences that users won't expect.

[jzemmels]

My intention here was to return something similar to what's returned by the other read_waterdata functions, hence the wrapping in sf::st_as_sf. Do we want the result to be a "pure" data.frame instead? Oh, you probably mean converting from a data.table. Nvm, I'll change that.

[ldecicco-USGS]

Oh, I just meant not a data.table, it can still be sf

[ldecicco-USGS]

When you get to this part, don't re-invent the wheel, just copy/paste the text from the API documentation.

[jzemmels]

I was leaving this in case we wanted to expand the get_params and related functions to pull from the /statistics Swagger docs. I've not scraped Swagger documentation before, so I'm not sure how easy that would be to generalize. I can just copy + paste the API docs for now.

[ldecicco-USGS]

The OGC service has a schema service that provides descriptions and whatnot via an API (I swooned 🤩 when I saw it). I wouldn't scrape Swagger though, more technical debt than what's saved.

[ldecicco-USGS]

In many arguments in dataRetrieval, the user can enter a number, character, vector, a date - and the code will convert it all to a character. So, you can check if it makes sense for a particular argument, but (a) don't let that code get ridiculously big trying to do everything and (b) see if the natural errors make sense as-is. If not, maybe try to clean up that error message. Usually (I don't know about here specifically), if a user puts in something really out of left field, the query will fail with a message that does imply their request was bad. That being said, the "waterdata" functions all have a similar input, at some point we can explore if a argument checker function could be written in a way that's more helpful.

[ldecicco-USGS]

Ideally we can figure out a way to not write these out by hand. On the off-chance they add some columns or change the function will start to fail. (and yes, similar hard-coding has caused some headaches in the past with dataRetrieval). We added some logic (I think just in the current PR that removes max_results) to move all "id" columns to the far right (because users don't want to see those big ol'hashes first) UNLESS they are special like monitoring_location_id.