Generate a multistar object

In this case, we are going to suppose that we have tables of facts and dimensions in tibble format, which we have imported into R. In particular, we have the tables mrs_fact_age, mrs_fact_cause, mrs_where, mrs_when and mrs_who, obtained from the example presented in package starschemar: Two fact tables that share two of the three dimensions.

library(tidyr)
library(geomultistar)

ms <- multistar() %>%
  add_facts(
    fact_name = "mrs_age",
    fact_table = mrs_fact_age,
    measures = "n_deaths",
    nrow_agg = "count"
  ) 

ms <- ms %>%
  add_facts(
    fact_name = "mrs_cause",
    fact_table = mrs_fact_cause,
    measures = c("pneumonia_and_influenza_deaths", "other_deaths"),
    nrow_agg = "nrow_agg"
  )

ms <- ms %>%
  add_dimension(
    dimension_name = "where",
    dimension_table = mrs_where,
    dimension_key = "where_pk",
    fact_name = "mrs_age",
    fact_key = "where_fk"
  )

ms <- ms %>%
  add_dimension(
    dimension_name = "when",
    dimension_table = mrs_when,
    dimension_key = "when_pk",
    fact_name = "mrs_age",
    fact_key = "when_fk",
    key_as_data = TRUE
  ) %>%
  add_dimension(
    dimension_name = "who",
    dimension_table = mrs_who,
    dimension_key = "who_pk",
    fact_name = "mrs_age",
    fact_key = "who_fk"
  )

ms <- ms %>%
  relate_dimension(dimension_name = "where",
                   fact_name = "mrs_cause",
                   fact_key = "where_fk") %>%
  relate_dimension(dimension_name = "when",
                   fact_name = "mrs_cause",
                   fact_key = "when_fk")

Define geographic dimensions and attributes

To define the dimensions and geographic attributes of a multistar object, we must define a geomultistar specialization from it, which allows to store this information.

We create a geomultistar object from a multistar one defining the names of the dimensions that contain geographic information. In the example only one dimension.

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

gms <-
  geomultistar(ms, geodimension = "where")

For each attribute of a geographic dimension that we want to use in queries, we can define a vector layer with which a relationship can be established using one or more attributes of the dimension.

gms <- gms %>%
  define_geoattribute(
    attribute = "city",
    from_layer = usa_cities,
    by = c("city" = "city", "state" = "state")
  ) 

For the city attribute, a relationship is defined with a vector layer in sf format (usa_cities), using the city and state attributes² that have the same name in the layer.

Sometimes there may be problems establishing the correspondence between the geographic attributes and the vector layer: Instances may remain unrelated. To detect these situations, we can query the rows that do not have associated geometry using the following function.

empty_city <- gms %>%
  get_empty_geoinstances(attribute = "city")

The result obtained is shown below.

In this case, for the unknown cities, their location has not been determined. There may be several because other geographic data of less granularity may be known.

In the same way, the relationship for county with the corresponding layer (usa_counties) is defined.

gms <- gms %>%
  define_geoattribute(
    attribute = "county",
    from_layer = usa_counties,
    by = c("county" = "county", "state" = "state")
  )  

We check if they have all been related.

empty_county <- gms %>%
  get_empty_geoinstances(attribute = "county")

And the result obtained is shown below.

It also happens for the same instances. In this case we can see that the associated geometry is of a different type.

In the case of state the definition is carried out by associating the code to the corresponding one in the layer (usa_states).

gms <- gms %>%
  define_geoattribute(
    attribute = c("state"),
    from_layer = usa_states,
    by = c("state" = "state")
  ) 

Additionally, for an attribute we can generate its layer from the one associated with another related attribute of the dimension. This is what has been done below for division.

gms <- gms %>%
  define_geoattribute(
    attribute = "division",
    from_attribute = "state"
  ) 

In this case, the vector layer is generated from the data available in the layer under consideration. Sometimes this is precisely what is desired. If not, look for a vector layer at that level of detail.

If no attribute name is indicated, this operation is applied to the rest of the attributes of the dimension that do not have an associated vector layer by any of the methods presented, as shown below.

gms <- gms %>%
  define_geoattribute(from_attribute = "state")

With this we have all the attributes of the dimension with an associated layer, defined at its level of granularity³. On the other hand, we can change the layer of any attribute at any time, independently of the rest.

Define multidimensional queries

Multidimensional queries are defined from a multistar object in the starschemar package. Since we are working with geomultistar, a specialization of multistar, we can use this object to define the query.

library(starschemar)

gdqr <- dimensional_query(gms) %>%
  select_dimension(name = "where",
                   attributes = c("division_name", "region_name")) %>%
  select_dimension(name = "when",
                   attributes = c("year", "week")) %>%
  select_fact(name = "mrs_age",
              measures = c("n_deaths")) %>%
  select_fact(
    name = "mrs_cause",
    measures = c("pneumonia_and_influenza_deaths", "other_deaths")
  ) %>%
  filter_dimension(name = "when", week <= "03")

The definition functions of this type of queries are quite intuitive, they are explained in the vignette of package starschemar.

If we use the function run_query of the mentioned package, a normal multidimensional query is executed, as shown below.

ft <- gdqr %>%
  run_query() %>%
  multistar_as_flat_table()

The result of this type of query is a multistar object, which we can transform into a flat table to display it more easily, as we do below for the first rows..

Run queries adding geographic information

If instead of executing the standard query, we execute run_geoquery function, from this package, we automatically obtain a vector layer at the finest possible level of detail, depending on the definition of the query.

vl_sf <- gdqr %>%
  run_geoquery()

The first rows of the result can be seen below in table form.

The result is a vector layer that we can save, perform spatial analysis or queries on it, or we can see it as a map, using the functions associated with the sf class.

class(vl_sf)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"

plot(vl_sf[,"n_deaths"])

Although we have indicated in the query the attributes division_name and region_name, as can be seen in the figure, the result obtained is at the finest granularity level, in this case at the division_name level.

Only the parts of the divisions made up of states where there is recorded data are shown. If we wanted to show the full extent of each division, we should have explicitly associated a layer at that level.

Get wide tables

In geographic layers, geographic objects usually are not repeated. The tables are wide: for each object the rest of the attributes are defined as columns. By means of the parameter wider we can indicate that we want a result of this type.

vl_sf_w <- gdqr %>%
  run_geoquery(wider = TRUE)

The first rows of the result can be seen below in table form.

We can see that the attributes that are multivalued for each geographic object have been eliminated from the result table, and new measurement columns have been generated: one for each combination of values of these attributes with the original measurements.

The meaning of the name of the columns of the measurements is part of the result obtained, also in table format, as can be seen below.

In this case there was only one variable with a multiplicity greater than one. If there were more variables in this situation, they would be added to this table in the same way.

This dictionary table and layer structure can be saved in GeoPackage format using the save_as_geopackage function.

save_as_geopackage(vl_sf_w, "division")

The GeoPackages thus obtained can be used directly, for example in QGIS.