1
\$\begingroup\$

This is a simple C-Squares encoder in R.

FYI, C-Squares is a geocoding algorithm that assigns GPS coordinates to a spatial grid with a chosen resolution. Its specifications is here: https://www.cmar.csiro.au/csquares/spec1-1.htm

My code only works with the resolution 0.05°×0.05° for the time being

# requires {dplyr}

# quadrant
Q_CSQ <- c("threshold" = 0L, "north east" = 1L, "south east" = 3L, "south west" = 5L, "north west" = 7L)
# sub-quadrant
sQ_CSQ <- c("threshold" = 5L, "south west" = 1L, "south east" = 2L, "north west" = 3L, "north east" = 4L)

#' @title identify quadrant
#' @param lon,lat GPS coords
#' @param quad 1 of the 2 quadrant vectors above
identify_quadrant <- function(lon, lat, quad) {
    threshold <- quad["threshold"]
    return(dplyr::case_when(
        lon <  threshold & lat >= threshold ~ quad["north west"],
        lon >= threshold & lat >= threshold ~ quad["north east"],
        lon >= threshold & lat <  threshold ~ quad["south east"],
        lon <  threshold & lat <  threshold ~ quad["south west"],
        TRUE                                ~ NA_integer_
    ))
}

#' @title extract tens, units, tenths, hundredths
#' @param x numeric vector
retrieve_decimal <- function(x) {
    remainder <- abs(x)

    tens <- as.integer(x / 10)
    remainder <- remainder - 10 * tens

    units <- as.integer(remainder)
    remainder <- remainder - units

    tenths <- as.integer(10 * remainder)
    remainder <- remainder - tenths / 10

    hundredths <- as.integer(100 * remainder)

    return(list("tens" = tens, "units" = units, "tenths" = tenths, "hundredths" = hundredths))
}

#' @title identify 3'×3' (0.05°×0.05°) C-SQUARES from GPS coordinates
#' @param LONGITUDE,LATITUDE GPS coords, accept vectors
#' @return cells id
identify_33cells_csquares <- function(LONGITUDE, LATITUDE) {

    # extract tens, units, tenths, hundredths
    vec_lon <- retrieve_decimal(LONGITUDE)
    vec_lat <- retrieve_decimal(LATITUDE)

    # quadrants
    quad_tens       <- identify_quadrant(LONGITUDE,          LATITUDE,           Q_CSQ)
    quad_units      <- identify_quadrant(vec_lon$units,      vec_lat$units,      sQ_CSQ)
    quad_tenths     <- identify_quadrant(vec_lon$tenths,     vec_lat$tenths,     sQ_CSQ)
    quad_hundredths <- identify_quadrant(vec_lon$hundredths, vec_lat$hundredths, sQ_CSQ)

    return(sprintf("%d%d%02d:%d%d%d:%d%d%d:%d",
        quad_tens,   vec_lat$tens,   vec_lon$tens,
        quad_units,  vec_lat$units,  vec_lon$units,
        quad_tenths, vec_lat$tenths, vec_lon$tenths,
        quad_hundredths
    ))
}

The C-Squares algorithm has some edge cases (see section 9.4 Exceptions in the specifications), but for my current projet data, I don't really need to treat them.

I'd like to add the ability to choose an other resolution, and possibly treating the edge cases.

Thank you.

\$\endgroup\$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.