Package {SUNGEO}

SUNGEO: Sub-National Geospatial Data Archive: Geoprocessing Toolkit

Description

Integrates spatially-misaligned GIS datasets across incompatible geographic units. Part of the Sub-National Geospatial Data Archive System. For the underlying methods, see Zhukov, Byers, Davidson, and Kollman (2024) "Integrating Data Across Misaligned Spatial Units," Political Analysis, Volume 32, Number 1, pp. 17-33 doi:10.1017/pan.2023.5.

Author(s)

Maintainer: Yuri M. Zhukov ymz2@georgetown.edu

Authors:

• Yuri M. Zhukov ymz2@georgetown.edu

• Jason Byers

• Marty Davidson

See Also

Useful links:

• https://www.sungeo.org/

• https://github.com/zhukovyuri/SUNGEO

SUNGEO

Description

Sub-National Geospatial Data Archive System: Geoprocessing Toolkit

Details

See the README on [GitHub](https://github.com/zhukovyuri/SUNGEO#readme)

Data availability through SUNGEO API

Description

Census of geospatial and processed data files available to download using SUNGEO::get_data().

Usage

data(available_data)

Format

List of 42 data.table objects Geoset:GADM :Classes ‘data.table’ and 'data.frame': 249 obs. of 4 variables Geoset:GAUL :Classes ‘data.table’ and 'data.frame': 242 obs. of 4 variables Geoset:geoBoundaries :Classes ‘data.table’ and 'data.frame': 197 obs. of 4 variables Geoset:GRED :Classes ‘data.table’ and 'data.frame': 74 obs. of 4 variables Geoset:HEXGRID :Classes ‘data.table’ and 'data.frame': 199 obs. of 4 variables Geoset:MPIDR :Classes ‘data.table’ and 'data.frame': 52 obs. of 4 variables Geoset:NHGIS :Classes ‘data.table’ and 'data.frame': 1 obs. of 4 variables Geoset:PRIOGRID :Classes ‘data.table’ and 'data.frame': 199 obs. of 4 variables Geoset:SHGIS :Classes ‘data.table’ and 'data.frame': 68 obs. of 4 variables

country_iso3

Codes for available countries (ISO 3166-1 alpha-3). Character string.

country_name

Names of available countries. Character string.

geoset_years

Years of available historical boundary files. Character string.

space_units

Available spatial units of analysis. Character string.

Elections:LowerHouse:CLEA :Classes ‘data.table’ and 'data.frame': 168 obs. of 6 variables Demographics:Ethnicity:EPR :Classes ‘data.table’ and 'data.frame': 180 obs. of 6 variables Demographics:Ethnicity:GREG :Classes ‘data.table’ and 'data.frame': 234 obs. of 6 variables Demographics:Population:GHS :Classes ‘data.table’ and 'data.frame': 257 obs. of 6 variables Events:PoliticalViolence:ABADarfur :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ACLED :Classes ‘data.table’ and 'data.frame': 100 obs. of 6 variables Events:PoliticalViolence:BeissingerProtest :Classes ‘data.table’ and 'data.frame': 15 obs. of 6 variables Events:PoliticalViolence:BeissingerRiot :Classes ‘data.table’ and 'data.frame': 15 obs. of 6 variables Events:PoliticalViolence:BeissingerUkraine :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:COCACW :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCAfghanistanWITS :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCIraqSIGACT :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCIraqWITS :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCMexicoDrugRelatedMurders :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCMexicoHomicide :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCPakistanBFRS :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:ESOCPakistanWITS :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:GED :Classes ‘data.table’ and 'data.frame': 121 obs. of 6 variables Events:PoliticalViolence:Lankina :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:NIRI :Classes ‘data.table’ and 'data.frame': 12 obs. of 6 variables Events:PoliticalViolence:NVMS :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:PITF :Classes ‘data.table’ and 'data.frame': 133 obs. of 6 variables Events:PoliticalViolence:SCAD :Classes ‘data.table’ and 'data.frame': 60 obs. of 6 variables Events:PoliticalViolence:yzCaucasus2000 :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:yzChechnya :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:yzLibya :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Events:PoliticalViolence:yzUkraine2014 :Classes ‘data.table’ and 'data.frame': 1 obs. of 6 variables Infrastructure:Roads:gRoads :Classes ‘data.table’ and 'data.frame': 240 obs. of 6 variables Infrastructure:NightLights:DMSP :Classes ‘data.table’ and 'data.frame': 257 obs. of 6 variables PublicHealth:Covid19:JHUCSSEC19 :Classes ‘data.table’ and 'data.frame': 207 obs. of 6 variables Terrain:Elevation:ETOPO1 :Classes ‘data.table’ and 'data.frame': 256 obs. of 6 variables Terrain:LandCover:GLCC :Classes ‘data.table’ and 'data.frame': 257 obs. of 6 variables Weather:AirTemperatureAndPrecipitation:NOAA :Classes ‘data.table’ and 'data.frame': 209 obs. of 6 variables

country_iso3

Codes for available countries (ISO 3166-1 alpha-3). Character string.

country_name

Names of available countries. Character string.

year_range

Range of available years for data topic. Character string.

time_units

Available time units. Character string.

space_units

Available spatial units. Character string.

geosets

Names of available geographic boundary data sources. Character string.

Source

Sub-National Geospatial Data Archive System: Geoprocessing Toolkit (updated March 17, 2023).

Country code dictionary

Description

Reference table of country names and ISO-3166 codes, adapted from countrycode package.

Usage

data(cc_dict)

Format

data.table object, with 8626 obs. of 3 variables:

country_name

Country names. Character string.

country_name_alt

Alternative spellings of country names, ASCII characters only. Character string.

country_iso3

Country codes (ISO 3166-1 alpha-3). Character string.

Source

Vincent Arel-Bundock. Package countrycode: Convert Country Names and Country Code, version 1.40. CRAN (October 12, 2022).

Constituency level results for lower chamber legislative elections, Germany 2009.

Description

A simple feature collection containing the spatial geometries of electoral constituency borders, and data on turnout levels, votes shares and other attributes of lower chamber legislative elections.

Usage

data(clea_deu2009)

Format

Simple feature collection with 16 features and 10 fields. geometry type: MULTIPOLYGON. dimension: XY. bbox: xmin: 5.867281 ymin: 47.27096 xmax: 15.04388 ymax: 55.05902. epsg (SRID): 4326. proj4string: +proj=longlat +datum=WGS84 +no_defs.

cst

Constituency number. Numeric.

cst_n

Constituency name. Character.

ctr

Country number. Numeric.

ctr_n

Country name. Character.

yrmo

Year and month of election (YYYYMM). Character.

to1

Turnout in first round. Numeric.

vv1

Number of valid votes in first round. Numeric.

pvs1_margin

Popular vote share margin in first round. Numeric.

incumb_pty_n

Incumbent party name.

win1_pty_n

Party name of popular vote share winner in first round. Character.

Source

Constituency-Level Elections Archive (CLEA) https://electiondataarchive.org/

Constituency level results for lower chamber legislative elections, Germany 2009.

Description

A data.frame object containing the geographic centroids of electoral contituencies, and data on turnout levels, votes shares and other attributes of lower chamber legislative elections.

Usage

data(clea_deu2009_df)

Format

data.frame with 16 observations and 12 variables.

cst

Constituency number. Numeric.

cst_n

Constituency name. Character.

ctr

Country number. Numeric.

ctr_n

Country name. Character.

yrmo

Year and month of election (YYYYMM). Character.

to1

Turnout in first round. Numeric.

vv1

Number of valid votes in first round. Numeric.

pvs1_margin

Popular vote share margin in first round. Numeric.

incumb_pty_n

Incumbent party name.

win1_pty_n

Party name of popular vote share winner in first round. Character.

longitude

Longitude of constituency centroid. Numeric.

latitude

Latitude of constituency centroid. Numeric.

Source

Constituency-Level Elections Archive (CLEA) https://electiondataarchive.org/

Constituency level results for lower chamber legislative elections, Germany 2009.

Description

A simple feature collection containing the geographic centroids of electoral contituencies, and data on turnout levels, votes shares and other attributes of lower chamber legislative elections.

Usage

data(clea_deu2009_pt)

Format

Simple feature collection with 16 features and 10 fields. geometry type: POINT. dimension: XY. bbox: xmin: 6.953882 ymin: 48.54535 xmax: 13.40315 ymax: 54.18635. epsg (SRID): 4326. proj4string: +proj=longlat +datum=WGS84 +no_defs.

cst

Constituency number. Numeric.

cst_n

Constituency name. Character.

ctr

Country number. Numeric.

ctr_n

Country name. Character.

yrmo

Year and month of election (YYYYMM). Character.

to1

Turnout in first round. Numeric.

vv1

Number of valid votes in first round. Numeric.

pvs1_margin

Popular vote share margin in first round. Numeric.

incumb_pty_n

Incumbent party name.

win1_pty_n

Party name of popular vote share winner in first round. Character.

Source

Constituency-Level Elections Archive (CLEA) https://electiondataarchive.org/

Convert data.frame object into simple features object

Description

Function takes in x-, y-coordinates, and a data.frame of variables (optional) and returns an SFC object

Usage

df2sf(
  x_coord,
  y_coord,
  input_data = NULL,
  file = NULL,
  n_max = Inf,
  start = 0,
  projection_input = "EPSG:4326",
  zero.policy = FALSE,
  show_removed = FALSE
)

Arguments

Value

If show_removed==FALSE, returns an sf object, with rows corresponding to non-usable coordinates removed. If show_removed==TRUE, returns a list, with an sf object (Spatial_Coordinates), and a vector of indices corresponding to non-usable coordinates removed (Removed_Rows).

Examples

# Coordinates supplied as vectors
data(clea_deu2009_df)
out_1 <- df2sf(x_coord=clea_deu2009_df$longitude,y_coord = clea_deu2009_df$latitude)
class(out_1)
plot(out_1$geometry)

# Coordinates supplied as column mames
out_2 <- df2sf(x_coord="longitude",y_coord ="latitude", input_data = clea_deu2009_df)
plot(out_2["geometry"])

# Load from external file

tmp <- tempfile()
write.csv(clea_deu2009_df,file=tmp)
out_3 <- df2sf(x_coord="longitude",y_coord ="latitude", file=tmp)
plot(out_3["geometry"])

Fix polygon geometries

Description

Function to check validity and fix broken geometries in simple features polygon objects

Usage

fix_geom(x, n_it = 10)

Arguments

Value

Returns a sf polygon object, with self-intersections and other geometry problems fixed.

Examples

# Fix geometries for a single dataset
data(clea_deu2009)
out_1 <- fix_geom(clea_deu2009)
all(sf::st_is_valid(out_1))

Geocode addresses and coordinates with OpenStreetMap

Description

Finds geographic coordinates of addresses and place names (forward geocoding), or converts longitude/latitude coordinates to place names and administrative units (reverse geocoding), using OpenStreetMap's Nominatim API.

Usage

geocode_osm(
  query,
  match_num = 1,
  return_all = FALSE,
  details = FALSE,
  reverse = FALSE,
  lon = NULL,
  lat = NULL,
  zoom = 10,
  user_agent = NULL
)

Arguments

Details

Note that the Nominatim Usage Policy stipulates an absolute maximum of 1 request per second (https://operations.osmfoundation.org/policies/nominatim/). For batch geocoding of multiple addresses, please use geocode_osm_batch.

Value

A data.frame object.

For forward geocoding (reverse=FALSE):

• query. User-supplied address query. Character string.

• osm_id. OpenStreetMap ID. Character string.

• address. OpenStreetMap address. Character string.

• longitude. Horizontal coordinate. Numeric.

• latitude. Vertical coordinate. Numeric.

If details=TRUE, also contains:

• osm_type. OpenStreetMap feature type. Character string.

• importance. Relevance of Nominatim match to query, from 0 (worst) to 1 (best). Numeric.

• bbox_ymin. Minimum vertical coordinate of bounding box. Numeric.

• bbox_ymax. Maximum vertical coordinate of bounding box. Numeric.

• bbox_xmin. Minimum horizontal coordinate of bounding box. Numeric.

• bbox_xmax. Maximum horizontal coordinate of bounding box. Numeric.

For reverse geocoding (reverse=TRUE):

• osm_name. Full display name from OpenStreetMap. Character string.

• osm_adm0_code. ISO country code. Character string.

• osm_adm0. Country name. Character string.

• osm_adm1. State or first-level administrative division. Character string.

• osm_adm2. District or county. Character string.

• osm_adm3. Municipality. Character string.

• osm_adm4. Borough, village or street. Character string.

Examples

# Geocode an address (top match only)
geocode_osm("Michigan Stadium")

# Return detailed results for top match
geocode_osm("Michigan Stadium", details = TRUE)

# Return detailed results for all matches
geocode_osm("Michigan Stadium", details = TRUE, return_all = TRUE)

# Reverse geocode a coordinate pair
geocode_osm(reverse = TRUE, lon = -83.74868, lat = 42.26587, zoom = 18)

Batch geocode addresses and coordinates with OpenStreetMap

Description

Finds geographic coordinates of multiple addresses and place names (forward geocoding), or converts multiple longitude/latitude coordinate pairs to place names and administrative units (reverse geocoding), using OpenStreetMap's Nominatim API.

Usage

geocode_osm_batch(
  query = NULL,
  delay = 1,
  return_all = FALSE,
  match_num = 1,
  details = FALSE,
  reverse = FALSE,
  lon = NULL,
  lat = NULL,
  zoom = 10,
  user_agent = NULL,
  verbose = FALSE
)

Arguments

Details

Wrapper function for geocode_osm. Because the Nominatim Usage Policy stipulates an absolute maximum of 1 request per second, this function facilitates batch geocoding by adding a small delay between queries (https://operations.osmfoundation.org/policies/nominatim/).

Value

A data.frame object.

For forward geocoding (reverse=FALSE):

• query. User-supplied address query. Character string.

• osm_id. OpenStreetMap ID. Character string.

• address. OpenStreetMap address. Character string.

• longitude. Horizontal coordinate. Numeric.

• latitude. Vertical coordinate. Numeric.

If details=TRUE, also contains:

• osm_type. OpenStreetMap feature type. Character string.

• importance. Relevance of Nominatim match to query, from 0 (worst) to 1 (best). Numeric.

• bbox_ymin. Minimum vertical coordinate of bounding box. Numeric.

• bbox_ymax. Maximum vertical coordinate of bounding box. Numeric.

• bbox_xmin. Minimum horizontal coordinate of bounding box. Numeric.

• bbox_xmax. Maximum horizontal coordinate of bounding box. Numeric.

For reverse geocoding (reverse=TRUE):

• osm_name. Full display name from OpenStreetMap. Character string.

• osm_adm0_code. ISO country code. Character string.

• osm_adm0. Country name. Character string.

• osm_adm1. State or first-level administrative division. Character string.

• osm_adm2. District or county. Character string.

• osm_adm3. Municipality. Character string.

• osm_adm4. Borough or village. Character string.

Examples

# Geocode multiple addresses (top matches only)
geocode_osm_batch(c("Ann Arbor", "East Lansing", "Columbus"))


# With progress reports
geocode_osm_batch(c("Ann Arbor", "East Lansing", "Columbus"), verbose = TRUE)


# Return detailed results for all matches
geocode_osm_batch(c("Ann Arbor", "East Lansing", "Columbus"),
  details = TRUE, return_all = TRUE)


# Reverse geocode multiple coordinates (city level)
geocode_osm_batch(
  reverse = TRUE,
  lon = c(-83.743, -84.556, -82.999),
  lat = c(42.278,  42.701,  39.961),
  zoom = 10
)

Download data from SUNGEO server

Description

Function to download data files through the SUNGEO API. Function produces a data.table object, corresponding to the user's choice of countries, topics, sources, and spatial and temporal units.

Usage

get_data(
  country_names = NULL,
  country_iso3 = NULL,
  geoset = "GADM",
  geoset_yr = 2018,
  space_unit = "adm1",
  time_unit = "year",
  topics = NULL,
  year_min = 1990,
  year_max = 2017,
  print_url = TRUE,
  print_time = TRUE,
  error_stop = FALSE,
  by_topic = TRUE,
  skip_missing = TRUE,
  cache_param = FALSE,
  short_message = TRUE
)

Arguments

Value

data.table object, with requested data from SUNGEO API.

See Also

get_info

Examples

# Single country, single topic
out_1 <- get_data(country_name="Afghanistan",topics="Demographics:Population:GHS")
out_1



out_2 <- get_data(
	country_name=c("Afghanistan","Moldova"),
	topics=c("Demographics:Ethnicity:EPR","Demographics:Population:GHS"))
out_2



# Other boundary sets, spatial and time units
out_3 <- get_data(
	country_name="Albania",
	topics="Weather:AirTemperatureAndPrecipitation:NOAA",
	geoset="GAUL",geoset_yr=1990,space_unit="adm2",time_unit="month",
	year_min=1990,year_max=1991)
out_3

Information on available SUNGEO data files

Description

This function reports the availability of data files on the SUNGEO server, searchable by country and topic.

Usage

get_info(country_names = NULL, country_iso3s = NULL, topics = NULL)

Arguments

Value

list object, with three slots: 'summary', 'topics', and 'geoset'.

See Also

get_data

Examples

# Get list of all available data
out_1 <- get_info()
out_1["summary"]
out_1["topics"]
out_1["geosets"]

# Get list of available data for a single country
out_2 <- get_info(country_names="Afghanistan")
out_2

# Get list of available data for a single topic
out_3 <- get_info(topics="Elections:LowerHouse:CLEA")
out_3

# Get list of available data for a multiple countries and topics
out_4 <- get_info(
                 country_names=c("Afghanistan","Zambia"),
                 topics=c("Elections:LowerHouse:CLEA","Events:PoliticalViolence:GED"))
out_4

Population count raster for Germany, 2010.

Description

2.5 arc-minute resolution raster of estimates of human population (number of persons per pixel), consistent with national censuses and population registers, for the year 2010.

Usage

data(gpw4_deu2010)

Format

class : SpatRaster dimensions : 186, 220, 1 (nrow, ncol, nlyr) resolution : 0.04166667, 0.04166667 (x, y) extent : 5.875, 15.04167, 47.29167, 55.04167 (xmin, xmax, ymin, ymax) coord. ref. : lon/lat WGS 84 (EPSG:4326) source(s) : memory name : gpw_v4_population_count_rev11_2010_2pt5_min min value : 0.00 max value : 92915.66

Source

Gridded Population of the World (GPW) v4: Population Count, v4.11 doi:10.7927/H4JW8BX5.

Hexagonal grid for Germany.

Description

Regular hexagonal grid of 0.5 degree diameter cells, covering territory of Germany (2020 borders).

Usage

data(hex_05_deu)

Format

Simple feature collection with 257 features and 3 fields. geometry type: POLYGON. dimension: XY. bbox: xmin: 5.375001 ymin: 46.76568 xmax: 15.375 ymax: 55.13726. epsg (SRID): 4326. proj4string: +proj=longlat +datum=WGS84 +no_defs.

HEX_ID

Unique cell identifier. Character.

HEX_X

Longitude of cell centroid. Numeric.

HEX_Y

Latitude of cell centroid. Numeric.

Source

SUNGEO

Roads polylines for Germany, 1992

Description

Roads thematic layer from Digital Chart of the World. Subset: divided multi-lane highways.

Usage

data(highways_deu1992)

Format

Simple feature collection with 1741 features and 5 fields. geometry type: MULTILINESTRING. dimension: XY. bbox: xmin: 5.750933 ymin: 47.58799 xmax: 14.75109 ymax: 54.80712 epsg (SRID): 4326. proj4string: +proj=longlat +datum=WGS84 +no_defs.

MED_DESCRI

Is the road a divided multi-lane highway with a median? Character string.

RTT_DESCRI

Primary or secondary route? Character string.

F_CODE_DES

Feature code description (road or trail). Character string.

ISO

ISO 3166-1 alpha-3 country code. Character string.

ISOCOUNTRY

Country name. Character string.

Source

Defense Mapping Agency (DMA), 1992. Digital Chart of the World. Defense Mapping Agency, Fairfax, Virginia. (Four CD-ROMs). Available through DIVA-GIS: https://diva-gis.org/data.html (accessed May 19, 2026).

Automatically calculate Local G hot spot intensity

Description

Function automatically calculates the Local G hot spot intensity measure for spatial points or spatial polygons. Uses RANN for efficient nearest neighbor calculation (spatial points only); users can specify the number of neighbors (k). Users can specify the neighborhood style (see spdep::nb2listw) with default being standardized weight matrix (W).

Usage

hot_spot(
  insert,
  variable = NULL,
  style = "W",
  k = 9,
  remove_missing = TRUE,
  NA_Value = 0,
  include_Moran = FALSE
)

Arguments

Value

If input is sf, SpatialPolygonsDataFrame or SpatialPointsDataFrame object, returns sf object with same geometries and columns as input, appended with additional column containing Local G estimates (LocalG). If input is RasterLayer object, returns RasterBrick object containing original values (Original) and Local G estimates (LocalG).

Examples

# Calculate Local G for sf point layer
data(clea_deu2009_pt)
out_1 <- hot_spot(insert=clea_deu2009_pt, variable = clea_deu2009_pt$to1)
class(out_1)
plot(out_1["LocalG"])

# Calculate Local G for sf polygon layer (variable as numeric vector)
data(clea_deu2009)
out_2 <- hot_spot(insert=clea_deu2009, variable = clea_deu2009$to1)
summary(out_2$LocalG)
plot(out_2["LocalG"])


# Calculate Local G for sf polygon layer (variable as column name)
out_3 <- hot_spot(insert=clea_deu2009, variable = "to1")
summary(out_3$LocalG)
plot(out_3["LocalG"])



# Calculate Local G for a SpatialPolygonsDataFrame object (variable as column name)
out_4 <- hot_spot(insert=as(clea_deu2009,"Spatial"), variable = "to1")
summary(out_4$LocalG)
plot(out_4["LocalG"])

Line-in-polygon analysis

Description

Function for basic geometry calculations on polyline features, within an overlapping destination polygon layer.

Usage

line2poly(
  linez,
  polyz,
  poly_id,
  measurez = c("length", "density", "distance"),
  outvar_name = "line",
  unitz = "km",
  reproject = TRUE,
  na_val = NA,
  verbose = TRUE
)

Arguments

Value

An sf polygon object, with summary statisics of linez features aggregated to the geometries of polyz.

If measurez = "lengths", contains fields with suffixes

• "_length". Sum of line lengths within each polygon, in km or other units supplied in unitz.

If measurez = "density", contains fields with suffixes

• "_length". Sum of line lengths within each polygon, in km or other units supplied in unitz.

• "_area". Area of each polygon, in km^2 or the square of linear units supplied in unitz.

• "_density". Sum of line lengths divided by area of each polygon, in km/km^2 or other units supplied in unitz.

If measurez = "distance", contains fields with suffixes

• "_distance". Distance from each polygon to nearest line feature, in km or other units supplied in unitz.

If measurez = c("length","density","distance") (default), contains all of the above.

Examples

# Road lengths, densities and distance from polygon to nearest highway
data(hex_05_deu)
data(highways_deu1992)
out_1 <- line2poly(linez = highways_deu1992,
                   polyz = hex_05_deu,
                   poly_id = "HEX_ID")
plot(out_1["line_length"])
plot(out_1["line_density"])
plot(out_1["line_distance"])


# Replace missing road lengths and densities with 0's, rename variables
out_2 <- line2poly(linez = highways_deu1992,
                   polyz = hex_05_deu,
                   poly_id = "HEX_ID",
                   outvar_name = "road",
                   na_val = 0)
plot(out_2["road_length"])
plot(out_2["road_density"])
plot(out_2["road_distance"])

Make date ticker

Description

Function to create a table of consecutive dates, in SUNGEO-compliant format.

Usage

make_ticker(
  date_min = 19000101,
  date_max = as.integer(gsub("-", "", as.Date(Sys.Date())))
)

Arguments

Value

data.table object, with seven columns:

• DATE. Date in YYYYMMDD format. Integer.

• DATE_ALT. Date in Date (YYYY-MM-DD) format. Date.

• TID. Date ID, in consecutive integer format. Integer.

• YRWK. Week in YYYYWW format. Integer.

• WID. Weed ID, in consecutive integer format. Integer.

• YRMO. Month in YYYYMM format. Integer.

• MID. Month ID, in consecutive integer format. Integer.

• YEAR. Year in YYYY format. Integer.

Examples

# All dates from January 1, 1900 to today
out_1 <- make_ticker()
out_1

# All dates from January 1, 1200 to today
out_2 <- make_ticker(date_min=12000101)
out_2

# All dates from January 1, 1500 to December 31, 1899
out_3 <- make_ticker(date_min=15000101, date_max=18991231)
out_3

Merge list of tables on common variable(s)

Description

Function that finds a set of common columns in a list of tables, and merges the tables on these columns.

Usage

merge_list(lst)

Arguments

Value

data.table object

Examples

# Merge list of three tables with different common variables
A <- data.table::data.table(month=month.name,year=rep(1991:1992,each=12),A=rnorm(24))
B <- data.table::data.table(year=c(1991,1992),B=rbeta(2,1,1))
C <- data.table::data.table(month=month.name,C=runif(12))

out_1 <- merge_list(list(A,B,C))
out_1

Relative scale and nesting coefficients

Description

Function to calculate relative scale and nesting metrics for changes of support from a source polygon layer to an overlapping (but spatially misaligned) destination polygon layer.

Usage

nesting(
  poly_from = NULL,
  poly_to = NULL,
  metrix = "all",
  tol_ = 0.001,
  by_unit = FALSE
)

Arguments

Details

Currently supported metrics (metrix) include:

• Relative scale ("rs"). Measures whether a change-of-support (CoS) task is one of aggregation or disaggregation, by calculating the share of source units that are smaller than destination units. Its range is from 0 to 1, where values of 1 indicate pure aggregation (all source units are smaller than destination units) and values of 0 indicate no aggregation (all source units are at least as large as destination units). Values between 0 and 1 indicate a hybrid (i.e. some source units are smaller, others are larger than target units).

• Relative nesting ("rn"). Measures how closely source and destination boundaries align, by calculating the share of source units that cannot be split across multiple destination units. Its range is from 0 to 1, where values of 0 indicate no nesting (every source unit can be split across multiple destination units) and values of 1 indicate full nesting (no source unit can be split across multiple destination units).

• Relative scale, symmetric ("rs_sym"). Alternative measure of "rs", which ranges from -1 to 1. It calculates a difference between two proportions: the share of source units that is smaller than destination units (i.e. "rs" from standpoint of source units), and the share that is larger (i.e. "rs" from standpoint of destination units). Values of -1 indicate pure disaggregation (all source units are larger than destination units), 1 indicates pure aggregation (all source units are smaller than destination units). Values of 0 indicate that all source units are the same size as target units.

• Relative nesting, symmetric ("rn_sym"). Alternative measure of "rn", which ranges from -1 to 1. It calculates a difference between two components: the nesting of source units within destination units (i.e. "rn" from standpoint of source units), and the nesting of destination units within source units (i.e. "rn" from standpoint of destination units. Values of 1 indicate that source units are perfectly nested within destination units; -1 indicates that destination units are perfectly nested within source units.

• Relative scale, alternative ("rs_alt"). Alternative measure of "rs", rescaled as a proportion of destination unit area. This measure can take any value on the real line, with positive values indicating aggregation and negative values indicating disaggregation.

• Relative nesting, alternative ("rn_alt"). Alternative measure of "rn", which places more weight on areas of maximum overlap. The main difference between this measure and "rn" is its use of the maximum intersection area for each source polygon instead of averaging over the quadratic term. Two sets of polygons are considered nested if one set is completely contained within another, with as few splits as possible. If none or only a sliver of a source polygon area falls outside a single destination polygon, those polygons are "more nested" than a case where half of a source polygon falls in destination polygon A and half falls into another polygon B.

• Relative scale, conditional ("rs_nn"). Alternative measure of "rs", calculated for the subset of source units that are not fully nested within destination units.

• Relative nesting, conditional ("rn_nn"). Alternative measure of "rn", calculated for the subset of source units that are not fully nested within destination units.

• Proportion intact ("p_intact"). A nesting metric that requires no area calculations at all. This measure ranges from 0 to 1, where 1 indicates full nesting (i.e. every source unit is intact/no splits), and 0 indicates no nesting (i.e. no source unit is intact/all are split).

• Proportion fully nested ("full_nest"). A stricter version of "p_intact". This measure ranges from 0 to 1, where 1 indicates full nesting (i.e. every source unit is intact/no splits AND falls completely inside the destination layer), and 0 indicates no nesting (i.e. no source unit is both intact and falls inside destination layer).

• Relative overlap ("ro"). Assesses extent of spatial overlap between source and destination polygons. This measure is scaled between -1 and 1. Values of 0 indicate perfect overlap (there is no part of source units that fall outside of destination units, and vice versa). Values between 0 and 1 indicate a "source underlap" (some parts of source polygons fall outside of destination polygons; more precisely, a larger part of source polygon area falls outside destination polygons than the other way around). Values between -1 and 0 indicate a "destination underlap" (some parts of destination polygons fall outside of source polygons; a larger part of destination polygon area falls outside source polygons than the other way around). Values of -1 and 1 indicate no overlap (all source units fall outside destination units, and vice versa). This is a theoretical limit only; the function returns an error if there is no overlap.

• Gibbs-Martin index of diversification ("gmi"). Inverse of "rn", where values of 1 indicate that every source unit is evenly split across multiple destination units, and 0 indicates that no source unit is split across any destination units.

It is possible to pass multiple arguments to metrix (e.g. metrix=c("rn","rs")). The default (metrix="all") returns all of the above metrics.

The function automatically reprojects source and destination geometries to Lambert Equal Area prior to calculation, with map units in meters.

Values of tol_ can be adjusted to increase or decrease the sensitivity of these metrics to small border misalignments. The default value discards polygon intersections smaller than 0.001 square meters in area.

Value

Named list, with numeric values for each requested metric in metrix. If by_unit==TRUE, last element of list is a data.table, with nesting metrics disaggregated by source unit, where the first column is a row index for the source polygon layer.

Examples

# Calculate all scale and nesting metrics for two sets of polygons
data(clea_deu2009)
data(hex_05_deu)
nest_1 <- nesting(
              poly_from = clea_deu2009,
              poly_to = hex_05_deu
              )
nest_1

# Calculate just Relative Nesting, in the opposite direction
nest_2 <- nesting(
              poly_from = hex_05_deu,
              poly_to = clea_deu2009,
              metrix = "rn"
              )
nest_2

Point-to-polygon interpolation, simple overlay method

Description

Function for assigning values from a source point layer to a destination polygon layer, using simple point-in-polygon overlays

Usage

point2poly_simp(
  pointz,
  polyz,
  varz,
  char_varz = NULL,
  funz = list(function(x) {
     sum(x, na.rm = TRUE)
 }),
  na_val = NA,
  drop_na_cols = FALSE
)

Arguments

Details

Assignment procedures are the same for numeric and character string variables. All variables supplied in varz are passed directly to the function specified in funz. If different sets of variables are to be aggregated with different functions, both varz and funz should be specified as lists (see examples below).

Value

Returns a sf polygon object, with variables from pointz assigned to the geometries of polyz.

Examples

# Assignment of a single variable (sums)
data(hex_05_deu)
data(clea_deu2009_pt)
out_1 <- point2poly_simp(pointz=clea_deu2009_pt,
                         polyz=hex_05_deu,
                         varz="vv1")
plot(out_1["vv1"])

# Replace NA's with 0's
out_2 <- point2poly_simp(pointz = clea_deu2009_pt,
                         polyz = hex_05_deu,
                         varz = "vv1",
                         na_val = 0)
plot(out_2["vv1"])


# Multiple variables, with different assignment functions
out_3 <- point2poly_simp(pointz = clea_deu2009_pt,
                         polyz = hex_05_deu,
                         varz = list(
                           c("to1","pvs1_margin"),
                           c("vv1"),
                           c("incumb_pty_n","win1_pty_n")),
                         funz = list(
                           function(x){mean(x,na.rm=TRUE)},
                           function(x){sum(x,na.rm=TRUE)},
                           function(x){paste0(unique(na.omit(x)),collapse=" | ") }),
                         na_val = list(NA_real_,0,NA_character_))
plot(out_3["pvs1_margin"])

Point-to-polygon interpolation, tessellation method

Description

Function for interpolating values from a source point layer to a destination polygon layer, using Voronoi tessellation and area/population weights.

Usage

point2poly_tess(
  pointz,
  polyz,
  poly_id,
  char_methodz = "aw",
  methodz = "aw",
  pop_raster = NULL,
  varz = NULL,
  pycno_varz = NULL,
  char_varz = NULL,
  char_assign = "biggest_overlap",
  funz = function(x, w) {
     stats::weighted.mean(x, w, na.rm = TRUE)
 },
  return_tess = FALSE,
  seed = 1
)

Arguments

Details

This function interpolates point data to polygons with a two-step process. In the first step (tessellation), each point is assigned a Voronoi cell, drawn such that (a) the distance from its borders to the focal point is less than or equal to the distance to any other point, and (b) no gaps between cells remain. The second step (interpolation) performs a polygon-in-polygon interpolation, using the Voronoi cells as source polygons.

Currently supported integration methods in the second step (methodz) include:

• Areal weighting ("aw"). Values from poly_from weighted in proportion to relative area of spatial overlap between source features and geometries of poly_to.

• Population weighting ("pw"). Values from poly_from weighted in proportion to relative population sizes in areas of spatial overlap between source features and geometries of poly_to. This routine uses a third layer (supplied in pop_raster) to calculate the weights.

When a list of variables are supplied and one methods argument specified, then the chosen method will be applied to all variables.

When a list of of variables are supplied and multiple methods arguments specified, then weighting methods will be applied in a pairwise order. For example, specifying varz = list(c("to1","pvs1_margin"), c("vv1")) and methodz = c('aw', 'pw') will apply areal weighting to the the first set of variables (to1 and pvs1_margin) and population weighing to the second set (vv1).

Interpolation procedures are handled somewhat differently for numeric and character string variables. For numeric variables supplied in varz, "aw" and/or "pw" weights are passed to the function specified in funz. If different sets of numeric variables are to be aggregated with different functions, both varz and funz should be specified as lists (see examples below).

For character string (and any other) variables supplied in char_varz, "aw" and/or "pw" weights are passed to the assignment rule(s) specified in char_assign. Note that the char_varz argument may include numerical variables, but varz cannot include character string variables.

Currently supported assignment rules for character strings (char_assign) include:

• "biggest_overlap". For each variable in char_varz, the features in poly_to are assigned a single value from overlapping poly_from features, corresponding to the intersection with largest area and/or population weight.

• "all_overlap". For each variable in char_varz, the features in poly_to are assigned all values from overlapping poly_from features, ranked by area and/or population weights (largest-to-smallest) of intersections.

It is possible to pass multiple arguments to char_assign (e.g. char_assign=c("biggest_overlap","all_overlap")), in which case the function will calculate both, and append the resulting columns to the output.

Value

If return_tess=FALSE, returns a sf polygon object, with variables from pointz interpolated to the geometries of polyz.

If return_tess=TRUE, returns a list, containing

• "result". The destination polygon layer. sf object.

• "tess". The (intermediate) Voronoi tessellation polygon layer. sf object.

Examples

# Interpolation of a single variable, with area weights
data(hex_05_deu)
data(clea_deu2009_pt)
out_1 <- point2poly_tess(pointz = clea_deu2009_pt,
                             polyz = hex_05_deu,
                             poly_id = "HEX_ID",
                             varz = "to1")
plot(out_1["to1_aw"])

# Extract and inspect tessellation polygons
out_2 <- point2poly_tess(pointz = clea_deu2009_pt,
                             polyz = hex_05_deu,
                             poly_id = "HEX_ID",
                             varz = "to1",
                             return_tess = TRUE)
plot(out_2$tess["to1"])
plot(out_2$result["to1_aw"])


# Interpolation of multiple variables, with area and population weights
data(gpw4_deu2010)
gpw4_deu2010 <- terra::rast(gpw4_deu2010) # unwrap PackedSpatRaster
out_3 <- point2poly_tess(pointz = clea_deu2009_pt,
                         polyz = hex_05_deu,
                         poly_id = "HEX_ID",
                         methodz = c("aw","pw"),
                         varz = list(
                           c("to1","pvs1_margin"),
                           c("vv1")
                         ),
                         pycno_varz = "vv1",
                         funz = list(
                           function(x,w){stats::weighted.mean(x,w)},
                           function(x,w){sum(x*w)}
                           ),
                         char_varz = c("incumb_pty_n","win1_pty_n"),
                         pop_raster = gpw4_deu2010)
plot(out_3["vv1_pw"])

Area and population weighted polygon-to-polygon interpolation

Description

Function for interpolating values from a source polygon layer to an overlapping (but spatially misaligned) destination polygon layer, using area and/or population weights.

Usage

poly2poly_ap(
  poly_from,
  poly_to,
  poly_to_id,
  geo_vor = NULL,
  methodz = "aw",
  char_methodz = "aw",
  pop_raster = NULL,
  varz = NULL,
  pycno_varz = NULL,
  char_varz = NULL,
  char_assign = "biggest_overlap",
  funz = function(x, w) {
     stats::weighted.mean(x, w, na.rm = TRUE)
 },
  seed = 1
)

Arguments

Details

Currently supported integration methods (methodz) include:

• Areal weighting ("aw"). Values from poly_from weighted in proportion to relative area of spatial overlap between source features and geometries of poly_to.

• Population weighting ("pw"). Values from poly_from weighted in proportion to relative population sizes in areas of spatial overlap between source features and geometries of poly_to. This routine uses a third layer (supplied in pop_raster) to calculate the weights.

It is possible to pass multiple arguments to methodz (e.g. methodz=c("aw","pw")), in which case the function will calculate both sets of weights, and append the resulting columns to the output.

Interpolation procedures are handled somewhat differently for numeric and character string variables. For numeric variables supplied in varz, "aw" and/or "pw" weights are passed to the function specified in funz. If different sets of numeric variables are to be aggregated with different functions, both varz and funz should be specified as lists (see examples below).

For character string (and any other) variables supplied in char_varz, "aw" and/or "pw" weights are passed to the assignment rule(s) specified in char_assign. Note that the char_varz argument may include numerical variables, but varz cannot include character string variables.

Currently supported assignment rules for character strings (char_assign) include:

• "biggest_overlap". For each variable in char_varz, the features in poly_to are assigned a single value from overlapping poly_from features, corresponding to the intersection with largest area and/or population weight.

• "all_overlap". For each variable in char_varz, the features in poly_to are assigned all values from overlapping poly_from features, ranked by area and/or population weights (largest-to-smallest) of intersections.

It is possible to pass multiple arguments to char_assign (e.g. char_assign=c("biggest_overlap","all_overlap")), in which case the function will calculate both, and append the resulting columns to the output.

Value

sf polygon object, with variables from poly_from interpolated to the geometries of poly_to.

Examples

# Interpolation of a single variable, with area weights
data(clea_deu2009)
data(hex_05_deu)
out_1 <- poly2poly_ap(poly_from = clea_deu2009,
              poly_to = hex_05_deu,
              poly_to_id = "HEX_ID",
              varz = "to1"
             )
plot(out_1["to1_aw"])


# Interpolation of multiple variables, with area weights
out_2 <- poly2poly_ap(
              poly_from = clea_deu2009,
              poly_to = hex_05_deu,
              poly_to_id = "HEX_ID",
              varz = list(
                c("to1","pvs1_margin"),
                c("vv1") ),
              pycno_varz = "vv1",
              funz = list(
                function(x,w){stats::weighted.mean(x,w)},
                function(x,w){sum(x*w)} ),
              char_varz = c("incumb_pty_n","win1_pty_n")
             )
plot(out_2["vv1_aw"])



# Interpolation of a single variable, with population weights
data(gpw4_deu2010)
gpw4_deu2010 <- terra::rast(gpw4_deu2010) # unwrap PackedSpatRaster
out_3 <- poly2poly_ap(poly_from = clea_deu2009,
                         poly_to = hex_05_deu,
                         poly_to_id = "HEX_ID",
                         varz = "to1",
                         methodz = "pw",
                         pop_raster = gpw4_deu2010)
plot(out_3["to1_pw"])



# Interpolation of a single variable, with area and population weights
out_4 <- poly2poly_ap(poly_from = clea_deu2009,
                         poly_to = hex_05_deu,
                         poly_to_id = "HEX_ID",
                         varz = "to1",
                         methodz = c("aw","pw"),
                         pop_raster = gpw4_deu2010)
plot(out_4["to1_aw"])

Convert simple features object into regularly spaced raster

Description

This function takes in an sf spatial object (polygon or point) and returns a regularly spaced RasterLayer. Reverse translation option allows users to create an sf polygon object from the regularly spaced RasterLayer. This function can also conver the sf object into a cartogram with a user-specified variable name.

Usage

sf2raster(
  polyz_from = NULL,
  pointz_from = NULL,
  input_variable = NULL,
  reverse = FALSE,
  poly_to = NULL,
  return_output = NULL,
  return_field = NULL,
  aggregate_function = list(function(x) mean(x, na.rm = TRUE)),
  reverse_function = list(function(x) mean(x, na.rm = TRUE)),
  grid_dim = c(1000, 1000),
  cartogram = FALSE,
  carto_var = NULL,
  message_out = TRUE,
  return_list = FALSE
)

Arguments

Value

If return_list=FALSE (default) and reverse=FALSE (default), returns RasterLayer object, with cell values corresponding to input_variable.

If return_list=TRUE and input layer is polygon, returns a list containing

• "return_output". Output raster, with values corresponding to input_variable. RasterLayer object.

• "return_centroid". Raster of centroids, with values corresponding to input_variable. RasterLayer object.

• "poly_to". Source polygons, with columns corresponding to input_variable and auto-generated numerical ID Field. sf object.

• "return_field". Output raster, with values corresponding to auto-generated numerical ID Field. RasterLayer object.

If return_list=TRUE and input layer is points, returns a list containing

• "return_output". Output raster, with values corresponding to input_variable. RasterLayer object.

• "return_point". Source points, with column corresponding to input_variable.

If reverse=TRUE, returns an sf polygon layer, with columns corresponding to input_variable and auto-generated numerical ID Field.

Examples

# Rasterization of polygon layer.
data(clea_deu2009)
out_1 <- sf2raster(polyz_from = utm_select(clea_deu2009),
                   input_variable = "to1")
terra::plot(out_1)


# Rasterization of point layer
data(clea_deu2009_pt)
out_2 <- sf2raster(pointz_from = utm_select(clea_deu2009_pt),
                   input_variable = "to1",
                   grid_dim = c(25,25))
terra::plot(out_2)


# Cartogram (vote turnout scaled by number of valid votes)
out_3 <- sf2raster(polyz_from = utm_select(clea_deu2009),
                   input_variable = "to1",
                   cartogram = TRUE,
                   carto_var = "vv1")
terra::plot(out_3)


# Polygonization of cartogram raster
out_4a <- sf2raster(polyz_from = utm_select(clea_deu2009),
                    input_variable = "to1",
                    cartogram = TRUE,
                    carto_var = "vv1",
                    return_list = TRUE)
out_4 <- sf2raster(reverse = TRUE,
                   poly_to = out_4a$poly_to,
                   return_output = out_4a$return_output,
                   return_field = out_4a$return_field)
terra::plot(out_4)

Smart numerical rounding function

Description

Function to round numerical values with minimal information loss (e.g. to avoid "0.000" values in tables).

Usage

smart_round(x, rnd = 0, return_char = TRUE)

Arguments

Details

Rounds the values in its first argument to the specified number of decimal places (default 0). If brute-force rounding produces zero values (e.g. "0.00"), the number of decimal places is expanded to include the first significant digit.

Value

If return_char=TRUE, returns a character string of same length as x. If return_char=FALSE, returns a numerical vector of same length as x.

Examples

# Round a vector of numbers, character string output (best for tables)
out_1 <- smart_round(c(.0013,2.3,-1,pi),rnd=2)
out_1

# Round a vector of numbers, numerical output
out_2 <- smart_round(c(.0013,2.3,-1,pi),rnd=2,return_char=FALSE)
out_2

Update bounding box of sf object

Description

Function to update the coordinates of the bounding box of sf vector data objects (e.g. after cropping or subsetting).

Usage

update_bbox(sfobj)

Arguments

Value

sf object, with corrected bounds.

Examples

# Update bbox for subset of sf object
data(clea_deu2009)
out_1 <- update_bbox(clea_deu2009[clea_deu2009$cst_n%in%c("Berlin"),])
out_1

# Bounding box of full dataset
data.table::as.data.table(clea_deu2009)[,sf::st_bbox(geometry)]

# Bounding box of subset (incorrect)
data.table::as.data.table(clea_deu2009)[cst_n%in%c("Berlin"),sf::st_bbox(geometry)]

# Corrected bounding box
data.table::as.data.table(out_1)[,sf::st_bbox(geometry)]

Automatically convert geographic (degree) to planar coordinates (meters)

Description

Function to automatically convert simple feature, spatial and raster objects with geographic coordinates (longitude, latitude / WGS 1984, EPSG:4326) to planar UTM coordinates. If the study region spans multiple UTM zones, defaults to Albers Equal Area.

Usage

utm_select(x, max_zones = 5, return_list = FALSE)

Arguments

Details

Optimal map projection for the object x is defined by matching its horizontal extent with that of the 60 UTM zones. If object spans multiple UTM zones, uses either the median zone (if number of zones is equal to or less than max_zones) or Albers Equal Area projection with median longitude as projection center (if number of zones is greater than max_zones).

Value

Re-projected layer. sf or RasterLayer object, depending on input.

If return_list=TRUE, returns a list object containing

• "x_out". The re-projected layer. sf or RasterLayer object, depending on input.

• "proj4_best".proj4string of the projection. Character string.

Examples

# Find a planar projection for an unprojected (WSG 1984) hexagonal grid of Germany
data(hex_05_deu)
sf::st_crs(hex_05_deu)
out_1 <- utm_select(hex_05_deu)
sf::st_crs(out_1)

# Find a planar projection for a raster
data(gpw4_deu2010)
gpw4_deu2010 <- terra::rast(gpw4_deu2010) # unwrap PackedSpatRaster
sf::st_crs(gpw4_deu2010)
out_2 <- utm_select(gpw4_deu2010)
sf::st_crs(out_2)