---
title: "Retail Location Analysis by Drive Time"
title-block: false
prefer-html: true

params:
  data_file: "Austin.xlsx"
  shopping_center_type: "Neighborhood Shopping Center" # Valid types = "Convenience Store", "Neighborhood Shopping Center", "Community Shopping Center", "Regional Shopping Center", and "Super-Regional Shopping Center"
  amenity_types:
    - school
    - college
    - university
    - restaurant
    - cafe
    - bank
    - pharmacy
    - supermarket
    - hospital
    - park
  weights:
    value:
      income: 0.30
      population_density: 0.30
      amenities: 0.10
      property_value: 0.10
      traffic_count: 0.20

format:
  html:
    toc: true
    toc-depth: 3
#    code-fold: true
#    code-tools: true
#    code-summary: "Show Code Block"
    theme: cosmo
    highlight-style: github
#  pdf:
#    toc: true
#    toc-depth: 3
#    number-sections: true
#    highlight-style: github
#    documentclass: article
#    geometry: margin=1in
#    mainfont: Calibri
#    pdf-engine: xelatex
#    header-includes:
#      - \usepackage{adjustbox}
#      - \usepackage{tabularx}
#      - \usepackage{booktabs}
#      - \usepackage{array}

execute:
  echo: false
  warning: false
  message: false
---
![](logo.png){width=260px style="float:left; margin-right:12px;"}

**Jim Amorin, CAE, MAI, SRA, AI-GRS, CDEI**  
Published: `r format(Sys.Date(), '%B %d, %Y')`

<br clear="all">

```{r}
#| include: false
options(
  kableExtra.latex.load_packages = FALSE,
  kableExtra.latex.table_env = "tabular",
  kableExtra.latex.longtable = FALSE,
  kableExtra.latex.bold_header = FALSE,
  knitr.kable.NA = ""
)
```

```{r}
#| label: load-libraries
#| include: false
#| message: false

# Core packages
library(tidyverse)        # Data manipulation and visualization
library(sf)               # Spatial data handling
library(tidycensus)       # Access to U.S. Census demographic data
library(janitor)          # Cleans and standardizes data frames
library(readxl)           # Reads Excel files (.xls and .xlsx)
library(writexl)          # Exports data frames to Excel format
library(scales)           # Formats numbers, percentages, and currency for display
library(httr)             # Makes HTTP requests for APIs (e.g., Census, Mapbox)
library(purrr)            # Functional programming tools for iteration and mapping
library(dplyr)            # Grammar of data manipulation
library(stringr)          # String operations and text processing
library(tibble)           # Enhanced data frame structure
library(plotly)           # Interactive visualizations and charts
library(osmdata)          # Retrieves spatial features from OpenStreetMap
library(purrr)            # Functional iteration (redundant; can remove if already loaded)

# Mapping packages
library(leaflet)          # Interactive web-based maps
library(leaflet.extras)   # Adds drawing tools, heatmaps, and search to leaflet maps
library(mapboxapi)        # Interfaces with Mapbox API for routing and isochrones

# Table formatting
library(knitr)            # Dynamic report generation for R Markdown / Quarto
library(kableExtra)       # Customizes and styles tables for HTML or PDF output



# Retrieve Mapbox token
mapbox_token <- Sys.getenv("MAPBOX_PUBLIC_TOKEN")

```

```{r}
#| label: geocoding
#| message: true
#| warning: false

library(httr)
library(dplyr)
library(readxl)
library(janitor)
library(stringr)

# ------------------------------------------------------------
# Load the Excel file and fill any missing coordinates in memory
# ------------------------------------------------------------

# 1. Read and clean the Excel file
properties <- read_excel(params$data_file) |> clean_names()

# 2. Identify columns
lat_col <- intersect(names(properties), c("lat", "latitude"))[1]
lon_col <- intersect(names(properties), c("lon", "longitude"))[1]
addr_fields <- intersect(names(properties), c("address", "city", "state", "zip", "zipcode"))

if (length(addr_fields) == 0)
  stop("No usable address fields found in the dataset.")

# Standardize names for simplicity
if (is.na(lat_col)) properties$lat <- NA_real_ else if (lat_col != "lat") properties <- rename(properties, lat = !!lat_col)
if (is.na(lon_col)) properties$lon <- NA_real_ else if (lon_col != "lon") properties <- rename(properties, lon = !!lon_col)

# 3. Define a simple geocoder
geocode_address <- function(address) {
  resp <- GET(
    "https://geocoding.geo.census.gov/geocoder/locations/onelineaddress",
    query = list(address = address, benchmark = "Public_AR_Current", format = "json"),
    timeout(15)
  )
  dat <- try(httr::content(resp, as = "parsed"), silent = TRUE)
  if (!inherits(dat, "try-error") && length(dat$result$addressMatches) > 0) {
    coords <- dat$result$addressMatches[[1]]$coordinates
    return(list(lat = coords$y, lon = coords$x))
  } else {
    return(list(lat = NA_real_, lon = NA_real_))
  }
}

# 4. Find and fill missing coordinates
missing_rows <- which(is.na(properties$lat) | is.na(properties$lon))
if (length(missing_rows) == 0) {
  message("✅ All coordinates present — no geocoding needed.")
} else {
  message("🌎 Found ", length(missing_rows), " record(s) missing coordinates. Geocoding now...")

  for (i in missing_rows) {
    full_address <- str_squish(paste(na.omit(unlist(properties[i, addr_fields])), collapse = ", "))
    coords <- geocode_address(full_address)

    # Explicitly write into the data frame and reassign
    properties$lat[i] <- coords$lat
    properties$lon[i] <- coords$lon

    message(sprintf("  • %s → lat %.6f, lon %.6f", full_address, coords$lat, coords$lon))
    Sys.sleep(1)
  }
}

# 5. Convert to spatial (for mapping later)
properties_sf <- properties |>
  filter(!is.na(lat) & !is.na(lon)) |>
  sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

# Keep both versions in memory
message("✅ properties (data frame) and properties_sf (sf object) ready in memory.")
```


# Introduction

This guide provides a comprehensive walk through of the Property Location Analysis Tool, an R-based system for evaluating real estate locations using geospatial data, demographics, and amenity analysis. This version uses drive-time isochrones instead of simple radius buffers to more accurately represent realistic trade areas based on actual road networks and travel times. Drive-time isochrones provide several advantages over circular radius buffers:

1.  Realistic trade areas: Isochrones follow actual road networks, accounting for highways, barriers, and traffic patterns
2.  Time-based accessibility: Reflects how far customers can actually travel in a given time period
3.  Irregular boundaries: Captures the natural asymmetry of accessible areas (e.g., faster travel along highways)
4.  Better demographic matching: More accurately identifies census tracts within the true catchment area

## Location Adjustments in Retail Property Appraisals

In the valuation of retail properties, location remains one of the most influential determinants of value. Location adjustments account for the qualitative and quantitative differences in the site characteristics that influence a property's utility, visibility, and access to a viable customer base. These adjustments are particularly critical when analyzing comparable sales in varying retail environments.

Retail success is closely tied to the property's trade area or that geographic region from which a retail establishment draws the majority of its customers. The size and shape of a trade area are rarely uniform and are influenced by multiple factors including existing and proposed residential developments, natural or manmade barriers (e.g., rivers, highways), and psychological boundaries such as perceptions of safety. Transportation infrastructure and commuter patterns further shape trade area dynamics by affecting ease of access for potential shoppers.

Demographic variables within a trade area such as population density, average household income, age distribution, and consumer spending patterns directly impact retail viability. For instance, household income levels and consumer expenditure habits (e.g., on groceries, apparel, or electronics) can indicate the buying power and preferences of the local market. Retailers often rely on these data points to determine store format, product mix, and even pricing strategies.

Retail property buyers and appraisers must also evaluate the visibility and accessibility of a site. Elements such as traffic volume, street frontage, line-of-sight from major roads, and access to public transportation all contribute to a location's competitive advantage. A site that is difficult to locate or reach whether due to excessive traffic speeds, poor signage, or other accessibility issues can significantly underperform in comparison to more accessible alternatives.

Ultimately, effective location adjustments require an understanding not only of physical proximity but of the socio-economic context within which a retail property operates. By carefully analyzing trade area boundaries using drive-time isochrones, demographic trends, and site-specific factors, appraisers can make more precise adjustments, thereby producing credible and supportable value conclusions.

@tbl-shopping-centers below outlines some of the key characteristics of various shopping center types, including their typical drive-time trade areas.

```{r}
#| label: tbl-shopping-centers
#| message: false
#| warning: false
#| echo: false

shopping_centers <- data.frame(
  Type = c(
    "Convenience Center",
    "Neighborhood Shopping Center",
    "Community Shopping Center",
    "Regional Shopping Center",
    "Super-Regional Shopping Center"
  ),
  Tenantry = c(
    "Stores that sell convenience goods (e.g., groceries, pharmaceuticals); not anchored by a supermarket.",
    "Stores that sell convenience goods and personal services (e.g., dry cleaning, shoe repair); often anchored by a supermarket.",
    "Stores that sell convenience goods, personal services, and shoppers' goods (e.g., apparel, appliances); may include a junior department store or off-price/discount store.",
    "Stores that sell general merchandise, shoppers' goods, and convenience goods; typically includes one or more department stores.",
    "Stores that sell general merchandise, apparel, furniture, home furnishings, services, and recreation; contains at least three major department stores."
  ),
  Size = c(
    "Less than 30,000 sq ft",
    "30,000 – 150,000 sq ft of gross leasable area; 4 – 10 acres",
    "100,000 – 300,000 sq ft of gross leasable area; 10 – 30 acres",
    "300,000 – 1,000,000 sq ft of gross leasable area; 30 – 100 acres",
    "Over 800,000 sq ft of gross leasable area"
  ),
  `Primary Trade Area` = c(
    "Less than 5-minute driving time",
    "Less than 5-minute driving time; 1 – 1½-mile range; 5,000 – 40,000 potential customers",
    "5 – 20-minute driving time; 3 – 6-mile range; 40,000 – 150,000 potential customers",
    "20 – 40-minute driving time; 5 – 10-mile range; 150,000 – 400,000 potential customers",
    "In excess of 30-minute driving time; typically 10 – 35-mile range; over 500,000 potential customers"
  ),
  stringsAsFactors = FALSE
)

# --- HTML version -----------------------------------------------------------
if (knitr::is_html_output()) {

  tbl <- kable(
    shopping_centers,
    caption = "Type of Shopping Centers",
    col.names = c("Type", "Tenantry", "Size", "Primary Trade Area"),
    booktabs = TRUE
  ) |>
    kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover"))

  note_html <- '<div style="text-align:left; font-size:85%; margin-top:4px;">
                  Source: Stephen F. Fanning, <i>Market Analysis for Real Estate</i>, 2nd ed. (Chicago: Appraisal Institute, 2014)
                </div>'
  knitr::asis_output(paste0(tbl, "\n", note_html))

# --- PDF (LaTeX) version ----------------------------------------------------
} else {

  latex_tbl <- "
\\begin{table}[htbp]
\\centering
\\caption{Type of Shopping Centers}
\\small
\\begin{tabular}{p{4cm} p{6cm} p{4cm} p{4cm}}
\\toprule
\\textbf{Type} & \\textbf{Tenantry} & \\textbf{Size} & \\textbf{Primary Trade Area} \\\\
\\midrule
Convenience Center & Stores that sell convenience goods (e.g., groceries, pharmaceuticals); not anchored by a supermarket. & Less than 30,000 sq ft & Less than 5-minute driving time \\\\
Neighborhood Shopping Center & Stores that sell convenience goods and personal services (e.g., dry cleaning, shoe repair); often anchored by a supermarket. & 30,000 – 150,000 sq ft of gross leasable area; 4 – 10 acres & Less than 5-minute driving time; 1 – 1½-mile range; 5,000 – 40,000 potential customers \\\\
Community Shopping Center & Stores that sell convenience goods, personal services, and shoppers' goods (e.g., apparel, appliances); may include a junior department store or off-price/discount store. & 100,000 – 300,000 sq ft of gross leasable area; 10 – 30 acres & 5 – 20-minute driving time; 3 – 6-mile range; 40,000 – 150,000 potential customers \\\\
Regional Shopping Center & Stores that sell general merchandise, shoppers' goods, and convenience goods; typically includes one or more department stores. & 300,000 – 1,000,000 sq ft of gross leasable area; 30 – 100 acres & 20 – 40-minute driving time; 5 – 10-mile range; 150,000 – 400,000 potential customers \\\\
Super-Regional Shopping Center & Stores that sell general merchandise, apparel, furniture, home furnishings, services, and recreation; contains at least three major department stores. & Over 800,000 sq ft of gross leasable area & In excess of 30-minute driving time; typically 10 – 35-mile range; over 500,000 potential customers \\\\
\\bottomrule
\\end{tabular}
\\vspace{2mm}
\\begin{flushleft}
\\footnotesize\\textit{Source: Stephen F. Fanning, Market Analysis for Real Estate, 2nd ed. (Chicago: Appraisal Institute, 2014)}
\\end{flushleft}
\\end{table}
"

  knitr::asis_output(latex_tbl)
}
```

## What This Tool Does

The Property Location Analysis Tool helps you determine a ranking of the location of comparables on a qualitative basis using a scoring of key elements:

-   **Provides an Input Section** to set parameters for the property type, weightings, etc.
-   **Score locations** based on multiple factors (income, population, amenities, property values, and traffic counts)
-   **Compare properties** side-by-side to identify the best opportunities
-   **Visualize results** with interactive maps, charts, and tables
-   **Generate reports** in multiple formats (Excel, HTML, PNG)
-   **Caches data** to avoid repeated API calls and speed up analysis

## Key Features

::: callout-note
## Core Capabilities

-   Geocoding addresses to coordinates (if necessary)
-   Fetching US Census demographic data by location
-   Querying OpenStreetMap for nearby amenities
-   Calculating weighted location scores
-   Creating interactive visualizations
:::

## Understanding the Code Structure

### Architecture Overview

The tool is organized into several functional modules as shown in the following graphic. :

::: {style="transform: scale(0.90); transform-origin: top center;"}
```{mermaid}
graph TD

A[Input: Property Details] --> B{Lat/Lon Available?}
B -- Yes --> C[Use Existing Coordinates]
B -- No --> D[Geocode Addresses]
C --> E[Data Collection]
D --> E
E --> F[Census Demographics]
E --> G[OSM Amenities]
F --> H[Scoring Engine]
G --> H
H --> I[Analysis Results]
I --> J[Visualizations]
I --> K[Reports]
J --> L[Maps, Charts, Tables]
K --> M[Excel, HTML]
```
:::

## Usage Instructions

### Basic Usage

1.  Prepare an Excel file with subject and sales' locations including these columns as a minimum:
    -   `Sale`: Identifies the sale number or Subject Property
    -   `Address`: Street address of the property
    -   `City`: City property lies in
    -   `State`: State location as two-letter abbreviation
    -   `Zip Code`: Zip Code for property location
    -   `Price/SF`: Price per square foot of property
    -   `Traffic Count`: For now you need to enter this figure from any available source
    -   Optional: `lat` or `latitude` column. **If you have it enter it** to avoid geocoding issues
    -   Optional: `lon` or `longitude` column. **If you have it enter it** to avoid geocoding issues
    -   Optional: other identifying columns
2.  Set parameters at the top of this Quarto document:
    -   `data_file`: Path to your Excel file
    -   `shopping_center_type`: Type of retail center
    -   `weights`: Set Custom Weights
3.  Render the document to generate analysis

# Setup and Configuration

### Mapbox API

To use drive-time isochrones, you need a **Mapbox Public Token key**:

1.  Create a free account at [mapbox.com](https://www.mapbox.com/)
2.  Generate an access token (refereed to in the code as MAPBOX_API_KEY)

```{r}
#| label: mapbox
#| eval: false


# Option 1: Set for current session
mb_access_token("MAPBOX_PUBLIC_TOKEN", install = FALSE)

# Option 2: Save to .Renviron (recommended)
# usethis::edit_r_environ()
```

### Census API Setup

```{r}
#| label: census-api
#| eval: false


# Get Census API key at: https://api.census.gov/data/key_signup.html
census_api_key("CENSUS_API_KEY", install = TRUE)
```

# Data Preparation

## Load Property Data

```{r}
#| label: load-data
#| message: false
#| warning: false
#| echo: false

# ------------------------------------------------------------
# Load and prepare dataset for downstream use
# ------------------------------------------------------------

# Confirm dataset exists
if (!exists("properties")) {
  stop("The 'properties' dataset was not found. Please run the geocoding chunk first.")
}

# Ensure required columns exist
required_cols <- c("sale", "lat", "lon")
missing_cols <- setdiff(required_cols, names(properties))
if (length(missing_cols) > 0) {
  stop("Missing expected columns in 'properties': ", paste(missing_cols, collapse = ", "))
}

# Keep ALL rows — even those missing coordinates
# (Do NOT filter out the Subject or others with missing lat/lon)
properties_all <- properties

# Create a spatial object using only those with valid coordinates
properties_sf <- properties_all |>
  dplyr::filter(!is.na(lat) & !is.na(lon)) |>
  sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

# Print confirmation
message("✅ Properties loaded: ", nrow(properties_all), " total (", nrow(properties_sf), " have coordinates).")

# Keep both data frames resident:
# - properties_all : full dataset, including Subject even if missing coordinates
# - properties_sf  : spatial object for mapping
```

```{r}
#| label: tbl-dataset
#| message: false
#| warning: false
#| echo: false

# ------------------------------------------------------------
# Display table of Subject and Comparable Sales (Geocoded)
# ------------------------------------------------------------

# Confirm dataset is available
if (!exists("properties")) {
  stop("The 'properties' dataset was not found. Please run the geocoding and load-data chunks first.")
}

# Select and display key columns
display_cols <- c("sale", "address", "city", "state", "lat", "lon")
display_cols <- intersect(display_cols, names(properties))

tbl_properties <- properties |>
  dplyr::select(all_of(display_cols)) |>
  dplyr::arrange(factor(sale, levels = unique(sale))) |>
  dplyr::mutate(
    lat = ifelse(is.na(lat), "— Missing —", round(as.numeric(lat), 6)),
    lon = ifelse(is.na(lon), "— Missing —", round(as.numeric(lon), 6))
  )

# Render table
knitr::kable(
  tbl_properties,
  caption = "Subject and Comparable Sales (Geocoded — Missing Coordinates Shown)",
  format = ifelse(knitr::is_latex_output(), "latex", "html"),
  align = c("l", "l", "l", "c", "r", "r")
)
```

# Drive Times Based on Shopping Center Type

```{r}
#| label: determine-drive-times


# Define drive time parameters based on shopping center type (single time per type)
drive_time_lookup <- tribble(
  ~type,                               ~primary_min,
  "Convenience Store",                 3,
  "Neighborhood Shopping Center",      5,
  "Community Shopping Center",         10,
  "Regional Shopping Center",          20,
  "Super-Regional Shopping Center",    30
)

# Get drive time for selected shopping center type
selected_drive_time <- drive_time_lookup |>
  filter(type == params$shopping_center_type) |>
  pull(primary_min)

if (length(selected_drive_time) == 0) {
  stop("Invalid shopping center type: ", params$shopping_center_type)
}

#cat("Shopping Center Type:", params$shopping_center_type, "\n")
#cat("Drive Time Isochrone:", selected_drive_time, "minutes\n")
```

This section determines drive time isochrones based on the the type of property under analysis. In this case, the property was defined as a `r params$shopping_center_type` which typically has a primary trade area of about `r selected_drive_time` minutes which will be used in this analysis. A review of @fig-map-isochrones displays these boundaries for each sale and the subject property.

```{r}
#| label: generate-isochrones
#| cache: true

sf_use_s2(FALSE)

# Rebuild spatial object to include any newly geocoded properties
properties_sf <- properties_all |>
  dplyr::filter(!is.na(lat) & !is.na(lon)) |>
  sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

# Ensure drive time is numeric
selected_drive_time <- as.numeric(selected_drive_time[1])

# Function to generate isochrones for a single property
get_isochrones <- function(lon, lat, time_min, profile = "driving") {
  tryCatch({
    message("Requesting isochrone for: ", round(lat, 6), ", ", round(lon, 6))
    iso <- mb_isochrone(
      location = c(as.numeric(lon), as.numeric(lat)),
      profile = profile,
      time = time_min
    )
    iso <- st_make_valid(iso)
    if (any(st_geometry_type(iso) == "GEOMETRYCOLLECTION"))
      iso <- st_collection_extract(iso, "POLYGON")
    iso |>
      mutate(
        origin_lon = lon,
        origin_lat = lat,
        time_minutes = time,
        area_sqkm = as.numeric(st_area(geometry)) / 1e6
      )
  }, error = function(e) {
    warning(paste("Failed to generate isochrone for", lon, lat, ":", e$message))
    return(NULL)
  })
}

# Generate isochrones for all properties
all_isochrones <- properties_sf |>
  mutate(
    lon = st_coordinates(geometry)[, 1],
    lat = st_coordinates(geometry)[, 2]
  ) |>
  st_drop_geometry() |>
  rowwise() |>
  mutate(isochrones = list(get_isochrones(lon, lat, selected_drive_time))) |>
  ungroup()

# Combine into one sf object
isochrones_sf <- all_isochrones |>
  filter(!sapply(isochrones, is.null)) |>
  unnest(isochrones) |>
  st_as_sf()
```

## Mapped Properties with Isochrones

```{r}
#| label: fig-map-isochrones
#| fig-cap: "Drive-Time Isochrones for Each Property"
#| warning: false
#| message: false
#| echo: false

# Identify subject and sales by the Sale column
subject_property <- properties_sf |> dplyr::filter(sale == "Subject")
sales_properties  <- properties_sf |> dplyr::filter(sale != "Subject")

# Label setup
subject_property <- subject_property |> dplyr::mutate(sale_label = "Subject")
sales_properties <- sales_properties |> 
  dplyr::mutate(sale_label = paste0("Sale ", sale))

# Combine and ensure sf geometry
properties_labeled <- dplyr::bind_rows(subject_property, sales_properties)
properties_labeled <- sf::st_as_sf(properties_labeled)

# Build map
map <- leaflet() |>
  addProviderTiles(providers$CartoDB.Positron) |>
  
  # Isochrones layer
  addPolygons(
    data = isochrones_sf,
    fillColor = "#D96C00",
    fillOpacity = 0.3,
    color = "#BF5700",
    weight = 2,
    label = ~paste0(sale, ": ", time_minutes, " min drive"),
    popup = ~paste0(
      "<b>", sale, "</b><br>",
      "Drive time: ", time_minutes, " minutes<br>",
      "Area: ", round(area_sqkm, 2), " sq km"
    ),
    group = "Isochrones",
    layerId = ~paste0("iso_", sale)  # Add layerId for tracking
  ) |>
  
  # Subject property (red)
  addCircleMarkers(
    data = subject_property,
    radius = 1.5,
    color = "darkred",
    fillColor = "red",
    fillOpacity = 0.9,
    group = "Subject Property",
    layerId = ~paste0("subj_", sale)
  ) |>
  
  # Comparable sales (blue)
  addCircleMarkers(
    data = sales_properties,
    radius = 1,
    color = "darkblue",
    fillColor = "blue",
    fillOpacity = 0.8,
    group = "Comparable Sales",
    layerId = ~paste0("sale_", sale)
  ) |>
  
  # Persistent text labels
  addLabelOnlyMarkers(
    data = properties_labeled,
    label = ~sale_label,
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(0, -10),
      style = list(
        "color" = "#000000",
        "font-weight" = "bold",
        "font-size" = "12px",
        "text-shadow" = "1px 1px 2px white"
      )
    ),
    group = "Labels"
  ) |>
  
  # Layer control
  addLayersControl(
    overlayGroups = c("Isochrones", "Subject Property", "Comparable Sales", "Labels"),
    options = layersControlOptions(collapsed = FALSE)
  ) |>
  
  # Add custom recenter button and enhanced functionality
  htmlwidgets::onRender("
    function(el, x) {
      var map = this;
      
      // ============================================
      // CREATE RECENTER BUTTON
      // ============================================
      var RecenterControl = L.Control.extend({
        options: {
          position: 'topleft'
        },
        onAdd: function(map) {
          var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control');
          var button = L.DomUtil.create('a', 'leaflet-control-recenter', container);
          button.innerHTML = '🎯';  // Bullseye emoji;
          button.href = '#';
          button.title = 'Recenter map to visible layers';
          button.style.fontSize = '20px';
          button.style.lineHeight = '30px';
          button.style.width = '30px';
          button.style.height = '30px';
          button.style.textAlign = 'center';
          button.style.textDecoration = 'none';
          button.style.backgroundColor = 'white';
          button.style.cursor = 'pointer';
          
          L.DomEvent.on(button, 'click', function(e) {
            L.DomEvent.preventDefault(e);
            recenterMap();
          });
          
          return container;
        }
      });
      
      map.addControl(new RecenterControl());
      
      // ============================================
      // RECENTER FUNCTION
      // ============================================
      function recenterMap() {
        var bounds = L.latLngBounds();
        var hasVisibleLayers = false;
        
        // Check each layer group
        map.eachLayer(function(layer) {
          // Check if layer is visible and has bounds
          if (map.hasLayer(layer) && layer.getBounds) {
            try {
              bounds.extend(layer.getBounds());
              hasVisibleLayers = true;
            } catch(e) {}
          }
          // Handle marker groups
          else if (map.hasLayer(layer) && layer.getLatLng) {
            try {
              bounds.extend(layer.getLatLng());
              hasVisibleLayers = true;
            } catch(e) {}
          }
        });
        
        // Fit bounds if we found visible layers
        if (hasVisibleLayers && bounds.isValid()) {
          map.fitBounds(bounds, {
            padding: [50, 50],
            maxZoom: 14
          });
        }
      }
      
      // ============================================
      // INITIAL FIT ON LOAD
      // ============================================
      setTimeout(function() {
        recenterMap();
      }, 100);
      
      // ============================================
      // RECENTER WHEN LAYERS TOGGLE
      // ============================================
      map.on('overlayadd overlayremove', function(e) {
        setTimeout(recenterMap, 50);
      });
      
      // ============================================
      // STYLING (includes your existing CSS)
      // ============================================
      var style = document.createElement('style');
      style.innerHTML = `
        /* Your existing layer control styles */
        .leaflet-control-layers-list label {
          display: flex;
          align-items: center;
          justify-content: flex-start;
          text-align: left !important;
          line-height: 1.2em;
          margin-bottom: 2px;
        }
        .leaflet-control-layers-list {
          padding-left: 6px;
          padding-right: 6px;
          font-size: 13px;
        }
        .leaflet-control-layers-selector {
          margin-right: 6px;
        }
        .leaflet-control-layers {
          background: rgba(255,255,255,0.95);
          box-shadow: 0 1px 3px rgba(0,0,0,0.3);
          border-radius: 6px;
        }
        
        /* Recenter button hover effect */
        .leaflet-control-recenter:hover {
          background-color: #f4f4f4;
        }
      `;
      el.appendChild(style);
    }
  ")

map
```

## Linking Isochrones and Census Data

The first major step in the analysis is to generate drive-time isochrones for each property. An isochrone defines the geographic area that can be reached within a specified travel time using the actual road network. This differs from a simple circular radius, which assumes straight-line distance and ignores highways, barriers, and travel speeds. By using drive-time boundaries, the analysis reflects a much more realistic picture of how far customers can travel to reach a property.

Once the isochrones are created, the next challenge is to associate those irregularly shaped areas with U.S. Census data. The Census Bureau provides demographic and economic data by fixed geographic units such as tracts, block groups, and blocks, but it does not allow direct data retrieval based on drive-time areas. In other words, while the Census API can return data within a given radius or within specific Census boundaries, it cannot natively return data for an isochrone polygon that follows real roads.

To bridge that gap, this tool takes a two-step approach. First, it downloads the full set of Census tracts for the state that contains the subject and comparable properties. Then, it uses a spatial intersection to determine which tracts fall wholly or partially within each drive-time isochrone. Because tract boundaries rarely align perfectly with isochrone shapes, the code calculates the proportion of each tract’s area that lies inside the isochrone and uses that percentage to weight demographic variables such as population and household income. This allows the results to reflect the demographic characteristics of the true trade area—based not on arbitrary circles but on how people actually travel.

In effect, the isochrone serves as a geographic “filter” that converts Census data from its static tract format into a dynamic representation of real-world accessibility. This integration of road-network travel times with demographic weighting provides a far more accurate foundation for evaluating retail trade areas and comparing locations on a consistent basis.

## Demographics Retrieval Function

A query to the US Census Bureau is necessary to fetch appropriate demographics for this analysis. As outlined, the key demographic variables used in this analysis include those shown in @tbl-census-vars.

```{r}
#| label: tbl-census-vars
#| tbl-cap: "Census Variables Used in Analysis"
#| echo: false

library(knitr)
library(kableExtra)

census_vars_df <- data.frame(
Variable = c("Population", "Median Income", "Median Home Value", "Median Gross Rent"),
Description = c(
"Total population in census tract",
"Median household income",
"Median value of owner-occupied homes",
"Median gross rent"
),
Census_Code = c("B01003_001", "B19013_001", "B25077_001", "B25046_001"),
stringsAsFactors = FALSE
)

# Rename the column for display

names(census_vars_df)[names(census_vars_df) == "Census_Code"] <- "Census Code"

kable(census_vars_df, align = c("l", "l", "c")) |>
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
```

```{r}
#| label: get-census-data
#| echo: false
#| message: false
#| warning: false
#| cache: true

# Get state from first property
state_abbr <- toupper(properties$state[1])

# Determine Census year
census_year <- 2023     # Latest year available

# Load census variables
# Updated to include Median Home Value and Median Gross Rent
census_vars <- c(
  total_population = "B01003_001",     # Total population
  median_income    = "B19013_001",     # Median household income
  median_home_value = "B25077_001",    # Median value of owner-occupied homes
  median_gross_rent = "B25046_001"     # Median gross rent
)

# Get census tracts for the state (suppress printed output)
suppressMessages({
  invisible(
    capture.output({
      tracts <- get_acs(
        geography = "tract",
        variables = census_vars,
        state = state_abbr,
        year = census_year,
        geometry = TRUE
      ) |>
        pivot_wider(
          id_cols = c(GEOID, NAME, geometry),
          names_from = variable,
          values_from = estimate
        ) |>
        st_as_sf()
    })
  )
})

# Transform to common CRS
tracts <- st_transform(tracts, st_crs(isochrones_sf))

# Make valid geometries
tracts <- st_make_valid(tracts)
# ------------------------------------------------------------------
# Export retrieved Census data to Excel
# ------------------------------------------------------------------
library(openxlsx)

# Define output filename (you can adjust path or naming)
output_file <- paste0("Census_Data_", state_abbr, "_", census_year, ".xlsx")

# Remove geometry before writing (Excel cannot store spatial objects)
tracts_table <- tracts |>
  st_drop_geometry()

# Write to Excel
write.xlsx(
  tracts_table,
  file = output_file,
  sheetName = "Census Data",
  overwrite = TRUE
)

cat("✅ Census data saved to", output_file, "\n")
```

## Spatial Join: Census Tracts Within Isochrones

This step identifies and summarizes demographic information within each property’s drive-time area (called an isochrone). Because Census data is published in tract boundaries rather than by drive-time zones, the analysis first determines which Census tracts fall inside each isochrone.

For each property, the code performs the following tasks:

1.  **Identify the Isochrone Area** The process begins by selecting the isochrone that corresponds to a specific property sale. Each isochrone represents the geographic area reachable by car within a set number of minutes from that location.

2.  **Find Census Tracts that Intersect the Isochron**e The program then looks for all Census tracts that overlap with the property’s isochrone boundary. Since these tracts rarely align perfectly with the drive-time area, only the overlapping portions are considered.

3.  **Weight Each Tract by Overlap Area** To ensure accuracy, the code calculates how much of each Census tract lies within the isochrone. It then weights the demographic data (population and income) based on the percentage of each tract included. For example, if only half of a tract falls inside the drive-time area, only half of that tract’s population and income values are counted.

4.  **Aggregate the Weighted Results** After weighting, the script adds up the adjusted population estimates and computes an average household income for the entire isochrone area. This produces a fair representation of the demographic characteristics actually within the drive-time zone rather than the full Census tracts.

5.  **Summarize by Property** Finally, the code compiles the results for each property, reporting the estimated total population, average income, and number of tracts that contributed data.

This analysis blends drive-time accessibility data with Census demographics. It provides localized, data-driven insight into the population and income patterns surrounding each property based on realistic travel areas rather than simple circular buffers.

```{r}
#| label: spatial-join
#| message: false
#| warning: false

get_census_in_isochrone <- function(sale_id) {

# Select the isochrone for the current sale

iso <- isochrones_sf |>
filter(sale == sale_id)

if (nrow(iso) == 0) {
return(tibble(
total_pop = NA_real_,
median_hh_income = NA_real_,
median_home_value = NA_real_,
median_gross_rent = NA_real_,
n_tracts = 0L
))
}

# Identify Census tracts intersecting the isochrone

tracts_in_iso <- st_intersection(tracts, iso) |>
mutate(
intersection_area = as.numeric(st_area(geometry)),
tract_area = as.numeric(st_area(st_geometry(tracts[match(GEOID, tracts$GEOID), ]))),
area_weight = intersection_area / tract_area,
weighted_pop = total_population * area_weight,
weighted_income = median_income * area_weight,
weighted_home_value = median_home_value * area_weight,
weighted_gross_rent = median_gross_rent * area_weight
)

# Aggregate weighted results

summary <- tracts_in_iso |>
st_drop_geometry() |>
summarise(
total_pop = sum(weighted_pop, na.rm = TRUE),
median_hh_income = weighted.mean(median_income, w = area_weight, na.rm = TRUE),
median_home_value = weighted.mean(median_home_value, w = area_weight, na.rm = TRUE),
median_gross_rent = weighted.mean(median_gross_rent, w = area_weight, na.rm = TRUE),
n_tracts = n()
)

return(summary)
}

# Compute Census summary for all properties

census_summary <- tibble(
sale = unique(isochrones_sf$sale)
) |>
rowwise() |>
mutate(census_data = list(get_census_in_isochrone(sale))) |>
unnest(census_data) |>
ungroup()
```

```{r}
#| label: tbl-census-summary
#| tbl-cap: "Demographic Summary within Drive-Time Isochrones"
#| echo: false
#| message: false
#| warning: false

library(knitr)
library(kableExtra)

if (exists("census_summary") && nrow(census_summary) > 0) {

census_summary_tbl <- census_summary |>
dplyr::select(
Sale = sale,
`Total Population` = total_pop,
`Median Household Income` = median_hh_income,
`Median Home Value` = median_home_value,
`Median Gross Rent` = median_gross_rent,
`Census Tracts Included` = n_tracts
) |>
dplyr::mutate(
`Total Population` = round(`Total Population`, 0),
`Median Household Income` = scales::dollar(`Median Household Income`, accuracy = 1),
`Median Home Value` = scales::dollar(`Median Home Value`, accuracy = 1),
`Median Gross Rent` = scales::dollar(`Median Gross Rent`, accuracy = 1)
)

kable(census_summary_tbl, align = c("l", "r", "r", "r", "r", "c")) |>
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
font_size = 12
)

} else {
cat("No census summary data available to display.\n")
}
```

## Calculate Population Density



```{r}
#| label: population-density
#| echo: false

# Add area (in square kilometers), convert to square miles, and calculate density

census_with_density <- census_summary |>
left_join(
isochrones_sf |>
st_drop_geometry() |>
select(sale, area_sqkm),
by = "sale"
) |>
mutate(
area_sqmi = area_sqkm * 0.386102,                # convert to square miles
pop_density_per_sqmi = total_pop / area_sqmi     # population per square mile
)
```

```{r}
#| label: tbl-pop-density
#| tbl-cap: "Population Density within Drive-Time Isochrones (per Square Mile)"
#| echo: false
#| message: false
#| warning: false

library(knitr)
library(kableExtra)
library(scales)

# Verify data exists

if (exists("census_with_density") && nrow(census_with_density) > 0) {

density_tbl <- census_with_density |>
mutate(
area_sqmi = round(area_sqmi, 2),
total_pop = round(total_pop, 0),
pop_density_per_sqmi = round(pop_density_per_sqmi, 1)
) |>
select(
Sale = sale,
`Total Population` = total_pop,
`Trade Area (sq mi)` = area_sqmi,
`Population Density (per sq mi)` = pop_density_per_sqmi
) |>
arrange(desc(`Population Density (per sq mi)`))

kable(
density_tbl,
align = c("l", "r", "r", "r"),
format = ifelse(knitr::is_latex_output(), "latex", "html")
) |>
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
font_size = 12
)

} else {
cat("No population density data available to display.\n")
}
```

@fig-pop-density-map visually reveals which trade areas contain the most concentrated residential populations. Areas shaded in darker tones correspond to higher population densities, indicating locations where retailers can expect stronger customer activity and greater potential sales volume within the same travel time. Lighter areas reflect lower residential concentration, signaling more dispersed customer bases or less intensive retail demand.

```{r}
#| label: fig-pop-density-map
#| fig-cap: "Population Density by Drive-Time Isochrone"
#| echo: false
#| message: false
#| warning: false

iso_with_density <- isochrones_sf |>
  left_join(
    census_with_density |>
      select(sale, pop_density_per_sqmi),
    by = "sale"
  )

subject_property <- properties_sf |> filter(sale == "Subject")
sales_properties <- properties_sf |> filter(sale != "Subject")

subject_property <- subject_property |> mutate(label = "Subject")
sales_properties <- sales_properties |> mutate(label = paste("Sale", sale))

# Improved color palette for better contrast on light basemaps
pal <- colorNumeric(
  palette = rev(viridisLite::plasma(8)),
  domain = iso_with_density$pop_density_per_sqmi
)

leaflet() |>
  addProviderTiles(providers$CartoDB.Positron) |>
  
  addPolygons(
    data = iso_with_density,
    fillColor = ~pal(pop_density_per_sqmi),
    color = "#444444",
    weight = 1,
    fillOpacity = 0.65,
    label = ~paste0(
      "Sale ", sale, "<br>",
      "Population Density: ", round(pop_density_per_sqmi, 1), " per sq mi"
    ),
    group = "Isochrones",
    layerId = ~paste0("iso_", sale)
  ) |>
  
  addCircleMarkers(
    data = subject_property,
    radius = 4,
    color = "darkred",
    fillColor = "red",
    fillOpacity = 0.9,
    group = "Subject",
    layerId = "subject_marker"
  ) |>
  
  addCircleMarkers(
    data = sales_properties,
    radius = 3,
    color = "black",
    fillColor = "black",
    fillOpacity = 0.8,
    group = "Sales",
    layerId = ~paste0("sale_", sale)
  ) |>
  
  addLabelOnlyMarkers(
    data = subject_property,
    label = ~label,
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(7, 0),
      style = list(
        "color" = "red",
        "font-weight" = "bold",
        "font-size" = "13px",
        "text-shadow" = "1px 1px 2px white"
      )
    )
  ) |>
  
  addLabelOnlyMarkers(
    data = sales_properties,
    label = ~label,
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(7, 0),
      style = list(
        "color" = "black",
        "font-weight" = "bold",
        "font-size" = "12px",
        "text-shadow" = "1px 1px 2px white"
      )
    )
  ) |>
  
  addLegend(
    pal = pal,
    values = iso_with_density$pop_density_per_sqmi,
    title = "Population Density<br>(per sq mi)",
    position = "bottomright"
  ) |>
  
  addLayersControl(
    overlayGroups = c("Isochrones", "Subject", "Sales"),
    options = layersControlOptions(collapsed = FALSE)
  ) |>
  
  # Add recenter button with automatic recentering
  htmlwidgets::onRender("
    function(el, x) {
      var map = this;
      
      // ============================================
      // CREATE RECENTER BUTTON
      // ============================================
      var RecenterControl = L.Control.extend({
        options: {
          position: 'topleft'
        },
        onAdd: function(map) {
          var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control');
          var button = L.DomUtil.create('a', 'leaflet-control-recenter', container);
          button.innerHTML = '🎯';  // Bullseye emoji;
          button.href = '#';
          button.title = 'Recenter map to visible layers';
          button.style.fontSize = '20px';
          button.style.lineHeight = '30px';
          button.style.width = '30px';
          button.style.height = '30px';
          button.style.textAlign = 'center';
          button.style.textDecoration = 'none';
          button.style.backgroundColor = 'white';
          button.style.cursor = 'pointer';
          
          L.DomEvent.on(button, 'click', function(e) {
            L.DomEvent.preventDefault(e);
            recenterMap();
          });
          
          return container;
        }
      });
      
      map.addControl(new RecenterControl());
      
      // ============================================
      // RECENTER FUNCTION
      // ============================================
      function recenterMap() {
        var bounds = L.latLngBounds();
        var hasVisibleLayers = false;
        
        map.eachLayer(function(layer) {
          if (map.hasLayer(layer) && layer.getBounds) {
            try {
              bounds.extend(layer.getBounds());
              hasVisibleLayers = true;
            } catch(e) {}
          }
          else if (map.hasLayer(layer) && layer.getLatLng) {
            try {
              bounds.extend(layer.getLatLng());
              hasVisibleLayers = true;
            } catch(e) {}
          }
        });
        
        if (hasVisibleLayers && bounds.isValid()) {
          map.fitBounds(bounds, {
            padding: [50, 50],
            maxZoom: 14
          });
        }
      }
      
      // ============================================
      // INITIAL FIT ON LOAD
      // ============================================
      setTimeout(function() {
        recenterMap();
      }, 100);
      
      // ============================================
      // RECENTER WHEN LAYERS TOGGLE
      // ============================================
      map.on('overlayadd overlayremove', function(e) {
        setTimeout(recenterMap, 50);
      });
      
      // ============================================
      // STYLING
      // ============================================
      var style = document.createElement('style');
      style.innerHTML = `
        .leaflet-control-recenter:hover {
          background-color: #f4f4f4;
        }
      `;
      el.appendChild(style);
    }
  ")
```

@fig-pop-density-bars highlights the relative concentration of residents within each property’s trade area. Higher bars correspond to locations with denser surrounding populations—an indicator of greater potential customer volume.”

```{r}
#| label: fig-pop-density-bars
#| fig-cap: "Relative Population Density by Property"
#| echo: false
#| message: false
#| warning: false

library(ggplot2)

ggplot(census_with_density, aes(x = reorder(sale, pop_density_per_sqmi), y = pop_density_per_sqmi)) +
geom_col(fill = "#BF5700") +
coord_flip() +
labs(
x = "Property",
y = "Population Density (per sq mi)",
title = "Relative Population Density by Trade Area"
) +
theme_minimal(base_size = 12)
```

# Amenity Analysis Using Isochrones

Queries OpenStreetMap [(https://wiki.openstreetmap.org/wiki/Key:amenity](https://wiki.openstreetmap.org/wiki/Key:amenity)) for nearby amenities within specified `r selected_drive_time` minute drive time for `r params$shopping_center_type`s. OSM’s `amenity` key is extremely broad. It includes everything from airports to post boxes.

In the context of retail demand and site selection, the relevant amenities are those that:

1.  Attract regular consumer visits (e.g., restaurants, cafes, banks, pharmacies).

2.  Represent daily or weekly necessities (e.g., supermarkets, convenience stores, schools).

3.  Encourage dwell time or foot traffic clustering (e.g., parks, places of worship, leisure).

4.  Support workforce presence (e.g., offices, post offices, government buildings).

These amenities serve as proxies for consumer presence, disposable income circulation, and accessibility to routine goods and services. For a retail valuation context, they also mirror the “market support” and “trade area vitality” concepts used in market analysis texts like Fanning, *Market Analysis for Real Estate*, 2nd ed.

The general amenity categories that are used in this analysis include the following:

-   Supermarkets - As a proxy for measuring the presence of food retail and daily household shopping activity.
-   Restaurants - Correlates to consumer engagement, social clustering and evening/weekend traffic.
-   Cafes - Also correlates to consumer engagement, social clustering and evening/weekend traffic.
-   Banks - Usually daytime traffic and draw to other nearby retail establishments.
-   Hospitals - Represents essential service and consistent daytime population.
-   Pharmacies - Regular essential service that often reflects high traffic.
-   School/College/University - Brings regular population flows including students, parents, and staff.
-   Parks - Encourages recreational foot traffic and neighborhood appeal.

```{r}
#| label: get-amenities
#| echo: false
#| cache: true

library(osmdata)

# Function to get amenities within isochrone
get_amenities_in_isochrone <- function(sale_id, amenity_type) {
  
  # Get isochrone boundary
  iso <- isochrones_sf |>
    filter(sale == sale_id)
  
  if (nrow(iso) == 0) {
    return(tibble(count = 0))
  }
  
  # Get bounding box
  bbox <- st_bbox(iso)
  
  # Query OSM
  amenities <- tryCatch({
    opq(bbox = bbox) |>
      add_osm_feature(key = "amenity", value = amenity_type) |>
      osmdata_sf()
  }, error = function(e) {
    return(NULL)
  })
  
  if (is.null(amenities)) {
    return(tibble(count = 0))
  }
  
  # Extract points
  amenity_points <- amenities$osm_points
  
  if (is.null(amenity_points) || nrow(amenity_points) == 0) {
    return(tibble(count = 0))
  }
  
  # Transform to match isochrone CRS
  amenity_points <- st_transform(amenity_points, st_crs(iso))
  
  # Filter to those within isochrone
  within_iso <- tryCatch({
    st_intersection(amenity_points, iso)
  }, error = function(e) {
    return(data.frame())
  })
  
  return(tibble(count = nrow(within_iso)))
}

# Get all amenities for all properties
amenity_counts <- expand_grid(
  sale = unique(isochrones_sf$sale),
  amenity_type = params$amenity_types
) |>
  rowwise() |>
  mutate(
    amenity_data = list(get_amenities_in_isochrone(sale, amenity_type))
  ) |>
  unnest(amenity_data) |>
  ungroup()

# Summarize total amenities per property
amenity_summary <- amenity_counts |>
  group_by(sale) |>
  summarise(
    total_amenities = sum(count, na.rm = TRUE),
    .groups = "drop"
  )

```

```{r}
#| label: tbl-amenities
#| tbl-cap: "Amenity Counts by Property"
#| echo: false
#| message: false
#| warning: false

library(dplyr)
library(knitr)
library(kableExtra)

# Pivot the amenity counts wider to show each amenity as a separate column

amenity_table <- amenity_counts |>
pivot_wider(
names_from = amenity_type,
values_from = count,
values_fill = 0
) |>
arrange(factor(sale, levels = unique(isochrones_sf$sale))) |>
mutate(
Total = rowSums(across(where(is.numeric)), na.rm = TRUE)
)

# Optional: order columns alphabetically except Total and Sale

cols_order <- c("sale", sort(setdiff(names(amenity_table), c("sale", "Total"))), "Total")
amenity_table <- amenity_table[, cols_order]

# Format table

kable(
amenity_table,
caption = "Amenity Counts by Property",
align = c("l", rep("r", ncol(amenity_table) - 1))
) |>
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
font_size = 12
)

```

# Traffic Count Integration

Traffic counts represent one of the most direct and quantifiable indicators of retail exposure. They measure the average number of vehicles or pedestrians passing a site over a given period, typically expressed as average daily traffic (ADT). For retailers, these counts serve as a proxy for visibility, accessibility, and potential customer flow.

High-traffic corridors generally translate to stronger impulse purchases, higher brand exposure, and greater potential for repeat visits. Conversely, properties with limited frontage on low-volume roads often face reduced visibility and a smaller trade area, even when nearby demographics appear favorable. Appraisers and market analysts use traffic data to evaluate a site’s competitive positioning, comparing comparable sales or lease rates within similar retail corridors to assess whether a location’s exposure justifies a premium or discount in value.

In short, traffic counts help quantify how often potential customers encounter a property in their daily travel patterns—bridging the gap between demographic potential and actual retail opportunity.

::: callout-note
Note: Traffic counts presently need to be added to the base Excel file as extraction of these variables is not presently a functionality.
:::

::: callout-note
## Traffic Count Data Sources

Traffic count data measure the average number of vehicles passing specific roadway points, often expressed as **Annual Average Daily Traffic (AADT)**. While this report uses data from the **Texas Department of Transportation (TxDOT)** [Open Data Portal](https://gis-txdot.opendata.arcgis.com/datasets/d5f56ecd2b274b4d8dc3c2d6fe067d37_0/explore), similar information is available from most U.S. states and many other countries.

### United States

-   **State Departments of Transportation (DOTs):**\
    Most state DOTs publish AADT datasets through their GIS or open data portals.\
    Examples include:
    -   [Caltrans (California)](https://data.ca.gov) → *Traffic Volumes*\
    -   [Florida Department of Transportation (FDOT)](https://geoplan.ufl.edu/) → *Traffic Counts*\
    -   [North Carolina DOT (NCDOT)](https://ncdot.maps.arcgis.com) → *Traffic Survey Data*
-   **Federal Highway Administration (FHWA):**\
    The FHWA provides nationwide roadway and traffic data through the\
    [Highway Performance Monitoring System (HPMS)](https://www.fhwa.dot.gov/policyinformation/hpms.cfm).

### Canada

-   **Provincial Ministries of Transportation:**\
    Each province maintains similar datasets:
    -   [Ontario Open Data Catalogue](https://data.ontario.ca) → *Annual Average Daily Traffic*\
    -   [DataBC Catalogue (British Columbia)](https://catalogue.data.gov.bc.ca) → *Traffic Counts*\
    -   [GeoDiscover Alberta](https://geodiscover.alberta.ca) → *Transportation / Traffic Volume Sites*
-   **Municipal and Regional Sources:**\
    Large municipalities (e.g., Toronto, Vancouver, Calgary) often publish local traffic counts through their open data portals.

### International

-   **United Kingdom:**\
    [Department for Transport Open Data](https://roadtraffic.dft.gov.uk/) — includes continuous and manual count sites across England, Scotland, and Wales.\
-   **European Union:**\
    Many member nations release traffic and transport statistics under\
    the [INSPIRE Directive](https://inspire.ec.europa.eu/) or through [Eurostat Transport Statistics](https://ec.europa.eu/eurostat/web/transport/data/database).\
-   **Global Alternatives:**\
    In regions without government data, commercial and open platforms such as\
    [TomTom Traffic](https://developer.tomtom.com/traffic-api/documentation/traffic-flow/traffic-flow-service),\
    [HERE Traffic](https://developer.here.com/documentation/traffic/dev_guide/index.html), or\
    [OpenStreetMap Traffic Layers](https://www.openstreetmap.org)\
    can provide approximate roadway congestion or flow estimates.

### *Future* Implementation Notes

To integrate local traffic data into this analysis: 1. Download the shapefile or GeoJSON version of your region’s AADT or traffic count dataset.\
2. Load it into R using `sf::st_read("path/to/file.shp")`.\
3. Perform a spatial join between your property points and the nearest road segment (`sf::st_join()` or `st_nearest_feature()`).\
4. Extract or average the relevant traffic count field for scoring.
:::

```{r}
#| label: tbl-traffic-counts
#| tbl-cap: "Traffic Counts by Property"
#| echo: false
#| message: false
#| warning: false

# Read the base Excel file defined in params

traffic_data <- read_excel(params$data_file) |> janitor::clean_names()

# Identify likely traffic count column (case-insensitive)

traffic_col <- names(traffic_data)[
str_detect(names(traffic_data), regex("traffic", ignore_case = TRUE))
][1]

if (is.na(traffic_col) || is.null(traffic_col)) {
stop("No column found in the Excel file that appears to contain traffic count data.")
}

# Try to find a sale identifier column (e.g., sale, id, site, property)

sale_col <- names(traffic_data)[
str_detect(names(traffic_data), regex("sale|id|site|property|subject", ignore_case = TRUE))
][1]

# If no sale column exists, create one using row numbers

if (is.na(sale_col) || is.null(sale_col)) {
traffic_data <- traffic_data |> mutate(sale = paste0("Site ", row_number()))
} else {
traffic_data <- traffic_data |> rename(sale = !!sale_col)
}

# Build the summary table

traffic_summary <- traffic_data |>
mutate(traffic_count = as.numeric(.data[[traffic_col]])) |>
select(sale, traffic_count)

# Compute summary metrics

avg_tc <- mean(traffic_summary$traffic_count, na.rm = TRUE)
max_tc <- max(traffic_summary$traffic_count, na.rm = TRUE)

# Append summary row (formatted text)

summary_row <- tibble(
sale = "Average / Max",
traffic_count = paste0(comma(round(avg_tc, 0)), " / ", comma(round(max_tc, 0)))
)

# Combine and format

traffic_summary_display <- traffic_summary |>
mutate(traffic_count = comma(traffic_count)) |>
bind_rows(summary_row)

# Rename columns for display

colnames(traffic_summary_display) <- c("Sale", "Traffic Count")

# Create formatted table

kable(
traffic_summary_display,
align = c("l", "r")
) |>
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
font_size = 12
)
```

## Weighting Function/Inputs

Required inputs to apply appropriate weight to the income (consumer spending power), the population density (a reflection of foot traffic potential), amenities (neighborhood quality and drawing power), property values (a function of costs to live in area), and traffic counts (a measure of visibility, accessibility and potential customers beyond population in trade area).

```{r}
#| label: tbl-weights
#| tbl-cap: "Model Weighting Parameters"
#| echo: false
#| message: false
#| warning: false

library(dplyr)
library(knitr)
library(kableExtra)
library(scales)

# Pull weighting values directly from YAML params

if (!is.null(params$weights$value)) {
weights <- unlist(params$weights$value)
} else {
warning("Weights not found in params — using default values.")
weights <- c(
income = 0.30,
population_density = 0.30,
amenities = 0.10,
property_value = 0.10,
traffic_count = 0.20
)
}

# Convert vector to a tidy data frame for display

weights_df <- tibble(
`Variable` = str_replace_all(names(weights), "_", " ") |> str_to_title(),
`Weight` = weights
) |>
mutate(
`Weight (%)` = percent(Weight, accuracy = 1)
) |>
select(`Variable`, `Weight (%)`)

# Create table

kable(
weights_df,
align = c("l", "r")
) |>
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
font_size = 12
)
```

## Scoring Function

Calculates a comprehensive location score (0-100) based on five factors:

::: panel-tabset
## Income Score

-   Formula: `rescale(median_income, to = c(0, 100), from = range(all_incomes))`
-   Higher income = better location (greater purchasing power)
-   Dyanmically normalized to 0-100 scale

**Interpretation:** The Income Score measures the relative purchasing power of households within each trade area. Median household income remains a reliable proxy for local spending capacity and the ability of residents to support a broad range of retail goods and services. Higher median incomes typically correspond with greater discretionary spending, stronger support for mid- and upper-tier retailers, and enhanced resilience during economic downturns.

This uses adaptive scaling to account for any range of income levels. The score is calculated by rescaling each trade area’s median income within the range of all incomes observed across the study sample. This approach preserves meaningful differences among higher-income areas, where a capped formula would otherwise treat them as identical, and allows for balanced comparison between affluent and moderate-income markets.

**Notes on the Scale:** The scaling is anchored to the lowest and highest median incomes within the comparison group. The location with the highest income receives a score of 100, while lowest income area scores a 0. Other income levels are scored proportionally.

## Population Density Score

-   Formula: `rescale(pop_density_per_sqmi, to = c(0, 100), from = range(all_densities))`
-   Higher density = more concentrated customer base and foot traffic potential
-   Normalized to 0-100 scale based on people per square mile

**Interpretation:** The Population Density Score represents the concentration of potential customers within the trade area, measured as population per square mile. Retail activity depends not only on total market size but also on proximity and frequency of visits. Areas with greater population density tend to generate stronger pass-by traffic, higher visibility, and more consistent sales volumes—particularly important for convenience-oriented retail like neighborhood shopping centers.

This metric uses adaptive scaling across the observed range of densities in the study sample. Properties in more densely populated areas score higher, reflecting their advantage in customer accessibility and repeat visit potential. A location with twice the density of another will have proportionally more scoring advantage, assuming other factors are equal.

**Notes on the Scale:** The lowest population density in the comparison group anchors the bottom of the scale at 0, while the highest density defines 100. All other sites are positioned proportionally between them based on their people-per-square-mile ratio.

## Amenities Score

-   Formula: `0.4 × (Diversity %) + 0.6 × (Density %, adaptively scaled 0–100)`
-   More diverse and concentrated amenities = stronger location
-   Balances variety and overall intensity of surrounding uses

**Interpretation:** The Amenities Score reflects both the diversity and density of nearby amenities that contribute to a site’s retail vitality. Diversity measures the proportion of tracked amenity categories—such as restaurants, schools, parks, banks, and pharmacies—that are represented within the trade area. Density measures the total number of amenities, scaled relative to the range observed across all study locations.

Together, these components create a balanced indicator of neighborhood convenience, consumer draw, and the underlying strength of the retail ecosystem. Locations that contain all tracked amenity types receive full credit for diversity, while those with fewer categories present earn proportionally less. The density component differentiates sites based on how intensively built and commercially active they are within the same amenity mix.

**Notes on the Scale:** The diversity index reaches 100 when all tracked amenity categories are represented within the trade area. The density index is adaptively normalized between the lowest and highest total amenity counts observed among all properties being compared, ensuring that each site’s score reflects its position within the actual range of market activity. The final score combines both effects, assigning 40% weight to diversity and 60% to density. This method ensures fair comparison across urban, suburban, and rural contexts without relying on arbitrary saturation thresholds.

## Property Value Score

-   Formula: `0.5 × rescale(median_home_value, 0-100) + 0.5 × rescale(median_gross_rent, 0-100)`
-   Higher values = more desirable neighborhood
-   Dynamically normalized to 0-100 scale based on observed range

**Interpretation:** The Property Value Score represents a composite measure of market vitality that combines indicators of long-term stability and short-term spending capacity within a trade area. Median home value reflects neighborhood affluence, ownership stability, and the general desirability of the surrounding residential market. Median gross rent complements this by capturing the relative cost of occupancy and the local balance between household income and housing expense—an important proxy for retail purchasing power and economic elasticity.

By weighting home value and rent equally (50/50), the score avoids bias toward either high-income owner-occupied areas or transient high-rent districts. A higher Property Value Score therefore signals a trade area that is both financially resilient and economically active, typically corresponding to locations that can support sustainable retail rents and consistent consumer demand.

High home value and high rent → strong, affluent market with high retail viability.\
High rent but low home value → possibly transient or overburdened market.\
High home value but low rent → stable ownership but low circulation spending.

**Notes on the Scale:** Both home values and rents are dynamically scaled to the 0-100 range based on the observed minimum and maximum within the comparison set. This ensures that the Property Value Score adapts to any market context—whether comparing properties in high-cost urban areas or more affordable suburban markets. The lowest home value in the dataset anchors the home value component at 0, while the highest anchors it at 100. The same logic applies independently to the rent component. This dynamic scaling maintains consistency with the Income Score and Population Density Score methodologies.

## Traffic Count Score

-   Formula: `rescale(traffic_count, to = c(0, 100), from = range(all_traffic_counts))`
-   Higher traffic = greater retail visibility and consumer flow
-   Dynamically normalized to 0–100 scale

**Interpretation:** The Traffic Count Score measures the retail exposure and accessibility of a location based on roadway volume. Higher traffic counts typically translate to greater visibility, impulse visits, and daily customer flow.

This score is dynamically scaled relative to the comparison group, ensuring that properties with exceptionally high or low counts remain distinguishable within each analysis set.

**Notes on the Scale:** The lowest observed traffic count anchors the scale at 0, and the highest at 100. When comparing sites across multiple markets, a shared range can be applied to maintain consistent interpretation across studies.
:::

**Total Score Calculation** $$
\text{Total Score} = \sum_{i=1}^{5} (\text{Component Score}_i \times \text{Weight}_i)
$$

**Interpretation:** The Total Score represents the weighted sum of five core indicators—Income, Population, Amenities, Property Value, and Traffic Count—each scaled from 0 to 100. Weights (𝑤ᵢ) control the relative importance of each factor and are defined in the document parameters. This additive framework allows the model to reflect both market potential (income, population), site quality (amenities, property value), and exposure/accessibility (traffic count).

**Notes on the Scale:** Each component score is dynamically normalized within its observed range to preserve proportional differences among locations. The composite score therefore expresses each site’s overall retail viability relative to all others in the analysis.

# Calculate Weighted Score for Location Analysis

Combine all metrics to generate location scores.

```{r}
#| label: calculate-scores
#| echo: false

# Combine all metrics (Census + Amenities)
location_metrics <- census_with_density |>
  left_join(amenity_summary, by = "sale") |>
  mutate(
    total_amenities = replace_na(total_amenities, 0),
    # Ensure density is available
    pop_density_per_sqmi = replace_na(pop_density_per_sqmi, 0)
  )

# Check for data
if (nrow(location_metrics) == 0) {
  stop("No location metrics data available.")
}

# Extract weights safely
if (!is.null(params$weights) && !is.null(params$weights$value)) {
  weights <- list(
    income = as.numeric(params$weights$value$income %||% 0.30),
    population_density = as.numeric(params$weights$value$population_density %||% 0.30),
    amenities = as.numeric(params$weights$value$amenities %||% 0.10),
    property_value = as.numeric(params$weights$value$property_value %||% 0.10),
    traffic_count = as.numeric(params$weights$value$traffic_count %||% 0.20)
  )
  message("✅ Using weights from YAML parameters.")
} else {
  warning("⚠️ Weights not found in params — using default values.")
  weights <- list(
    income = 0.30,
    population_density = 0.30,
    amenities = 0.10,
    property_value = 0.10,
    traffic_count = 0.20
  )
}

# Calculate Amenity Diversity (% of amenity types present)
amenity_diversity <- amenity_counts |>
  group_by(sale) |>
  summarise(
    unique_types = n_distinct(amenity_type[count > 0]),
    total_types = n_distinct(amenity_type),
    diversity_pct = 100 * unique_types / total_types,
    .groups = "drop"
  )

# Join diversity into metrics
location_metrics <- location_metrics |>
  left_join(amenity_diversity, by = "sale") |>
  mutate(diversity_pct = replace_na(diversity_pct, 0))

# Compute scores
location_scores <- location_metrics |>
  mutate(
    # 1. Income Score (Median Household Income)
    income_score = rescale(
      median_hh_income,
      to = c(0, 100),
      from = range(median_hh_income, na.rm = TRUE)
    ),

    # 2. Population Density Score (per square mile)
    density_score = rescale(
      pop_density_per_sqmi,
      to = c(0, 100),
      from = range(pop_density_per_sqmi, na.rm = TRUE)
    ),

    # 3. Amenity Score = 0.4 * diversity + 0.6 * amenity count
    amenity_density_scaled = rescale(
      total_amenities,
      to = c(0, 100),
      from = range(total_amenities, na.rm = TRUE)
    ),
    amenity_score = 0.4 * diversity_pct + 0.6 * amenity_density_scaled,

    # 4. Property Value Score (dynamically scaled)
    home_value_score = rescale(
      median_home_value,
      to = c(0, 100),
      from = range(median_home_value, na.rm = TRUE)
    ),
    rent_score = rescale(
      median_gross_rent,
      to = c(0, 100),
      from = range(median_gross_rent, na.rm = TRUE)
    ),
    property_value_score = 0.5 * home_value_score + 0.5 * rent_score,

    # 5. Traffic Count Score (properly matched by sale)
    traffic_score = {
      if ("traffic_count" %in% names(properties)) {
        # Create lookup and join to ensure proper alignment
        traffic_lookup <- properties |>
          select(sale, traffic_count)
        
        current_with_traffic <- location_metrics |>
          left_join(traffic_lookup, by = "sale")
        
        # Rescale based on matched data
        rescale(
          current_with_traffic$traffic_count,
          to = c(0, 100),
          from = range(current_with_traffic$traffic_count, na.rm = TRUE)
        )
      } else {
        rep(50, nrow(location_metrics))
      }
    }
  ) |>
  mutate(
    # 6. Weighted total score
    total_score = 
      income_score * weights$income +
      density_score * weights$population_density +
      amenity_score * weights$amenities +
      property_value_score * weights$property_value +
      traffic_score * weights$traffic_count
  ) |>
  arrange(desc(total_score)) |>
  mutate(rank = row_number())
```

```{r}
#| label: tbl-ranked-locations
#| tbl-cap: "Ranked Retail Locations"
#| echo: false
#| warning: false
#| message: false

# Confirm that location_scores exists

if (!exists("location_scores")) {
stop("Object 'location_scores' not found — run the calculate-scores chunk first.")
}

# Create ranked summary table

ranked_tbl <- location_scores |>
arrange(desc(total_score)) |>
mutate(Rank = row_number()) |>
select(
Rank,
Sale = sale,
`Total Score` = total_score,
`Income Score` = income_score,
`Population Score` = density_score,
`Amenities Score` = amenity_score,
`Property Value Score` = property_value_score,
`Traffic Score` = traffic_score
)

# Format and display

kable(
ranked_tbl,
digits = 1,
align = c("c", "l", rep("r", 6))
) |>
kable_styling(
full_width = TRUE,
position = "center",
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
font_size = 12
)
```

```{r}
#| label: attach-coordinates
#| echo: false
#| message: false
#| warning: false

# Reattach coordinates to location_scores

location_scores <- location_scores |>
left_join(
properties |>
select(sale, lat, lon),
by = "sale"
)
```

# Visualization: Scored Properties with Isochrones

```{r}
#| label: fig-final-map
#| fig-cap: "Property Locations Colored by Score with Drive-Time Isochrones"

# Color palette for scores
score_pal <- colorNumeric(
  palette = "RdYlGn",
  domain = location_scores$total_score
)

# Get top property isochrones
top_property <- location_scores$sale[1]
top_isochrones <- isochrones_sf |>
  filter(sale == top_property)

leaflet() |>
  addProviderTiles(providers$CartoDB.Positron) |>
  
  # Add isochrones for all properties (subtle)
  addPolygons(
    data = isochrones_sf,
    fillColor = "lightblue",
    fillOpacity = 0.1,
    color = "gray",
    weight = 1,
    group = "All Isochrones",
    layerId = ~paste0("iso_all_", sale)
  ) |>
  
  # Highlight top property isochrones
  addPolygons(
    data = top_isochrones,
    fillColor = "darkorange",
    fillOpacity = 0.3,
    color = "darkorange",
    weight = 2,
    label = ~paste0(time_minutes, " min drive"),
    group = "Top Property",
    layerId = ~paste0("iso_top_", sale)
  ) |>
  
  # Add scored properties
  addCircleMarkers(
    data = location_scores |>
      st_as_sf(coords = c("lon", "lat"), crs = 4326),
    radius = 8,
    fillColor = ~score_pal(total_score),
    fillOpacity = 0.8,
    color = "black",
    weight = 2,
    label = ~paste0(
      sale, " (Rank #", rank, "): ",
      round(total_score, 1), " points"
    ),
    popup = ~paste0(
      "<b>", sale, "</b><br>",
      "Rank: ", rank, "<br>",
      "Total Score: ", round(total_score, 1), "<br>",
      "Income Score: ", round(income_score, 1), "<br>",
      "Density Score: ", round(density_score, 1), "<br>",
      "Amenity Score: ", round(amenity_score, 1)
    ),
    layerId = ~paste0("score_", sale)
  ) |>
  
  addLegend(
    pal = score_pal,
    values = location_scores$total_score,
    title = "Location Score",
    position = "bottomright"
  ) |>
  
  addLayersControl(
    overlayGroups = c("All Isochrones", "Top Property"),
    options = layersControlOptions(collapsed = FALSE)
  ) |>
  
  # Add recenter button with automatic recentering
  htmlwidgets::onRender("
    function(el, x) {
      var map = this;
      
      // ============================================
      // CREATE RECENTER BUTTON
      // ============================================
      var RecenterControl = L.Control.extend({
        options: {
          position: 'topleft'
        },
        onAdd: function(map) {
          var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control');
          var button = L.DomUtil.create('a', 'leaflet-control-recenter', container);
          button.innerHTML = '🎯';  // Bullseye emoji;
          button.href = '#';
          button.title = 'Recenter map to visible layers';
          button.style.fontSize = '20px';
          button.style.lineHeight = '30px';
          button.style.width = '30px';
          button.style.height = '30px';
          button.style.textAlign = 'center';
          button.style.textDecoration = 'none';
          button.style.backgroundColor = 'white';
          button.style.cursor = 'pointer';
          
          L.DomEvent.on(button, 'click', function(e) {
            L.DomEvent.preventDefault(e);
            recenterMap();
          });
          
          return container;
        }
      });
      
      map.addControl(new RecenterControl());
      
      // ============================================
      // RECENTER FUNCTION
      // ============================================
      function recenterMap() {
        var bounds = L.latLngBounds();
        var hasVisibleLayers = false;
        
        map.eachLayer(function(layer) {
          if (map.hasLayer(layer) && layer.getBounds) {
            try {
              bounds.extend(layer.getBounds());
              hasVisibleLayers = true;
            } catch(e) {}
          }
          else if (map.hasLayer(layer) && layer.getLatLng) {
            try {
              bounds.extend(layer.getLatLng());
              hasVisibleLayers = true;
            } catch(e) {}
          }
        });
        
        if (hasVisibleLayers && bounds.isValid()) {
          map.fitBounds(bounds, {
            padding: [50, 50],
            maxZoom: 14
          });
        }
      }
      
      // ============================================
      // INITIAL FIT ON LOAD
      // ============================================
      setTimeout(function() {
        recenterMap();
      }, 100);
      
      // ============================================
      // RECENTER WHEN LAYERS TOGGLE
      // ============================================
      map.on('overlayadd overlayremove', function(e) {
        setTimeout(recenterMap, 50);
      });
      
      // ============================================
      // STYLING
      // ============================================
      var style = document.createElement('style');
      style.innerHTML = `
        .leaflet-control-recenter:hover {
          background-color: #f4f4f4;
        }
      `;
      el.appendChild(style);
    }
  ")
```

## Sensivity Analysis

The sensitivity analysis tests how stable each property’s ranking remains when the weighting of input factors such as income, population density, amenities, property values, and traffic counts is allowed to vary randomly within reasonable limits. This scenario analysis runs 5,000 times with alternative weighting configuration being simulated. Each configuration produces a new set of location rankings. The resulting distribution of simulated ranks for each property reflects how dependent its position is on the specific weighting assumptions used in the model.

The median simulated rank represents the property’s typical performance across all simulations, while the 5th and 95th percentiles mark its best-case and worst-case outcomes, respectively. A property with a low median rank and a narrow percentile range demonstrates consistent strength and resilience to changing assumptions. In contrast, a wide percentile range indicates that the property’s relative standing shifts considerably under different weighting scenarios, suggesting greater sensitivity to how the model prioritizes its components. Together, these measures reveal not only which locations perform best on average but also which are most robust to uncertainty in the underlying valuation logic.

::: callout-note
### About the Sensitivity Simulation

The “reasonable limits” for weighting variation are defined by a distribution, which produces random, positive weights that sum to one. This distribution allows each factor including Income, Population, Amenities, Property Value, and Traffic to fluctuate broadly between near-zero and moderate dominance, representing plausible shifts in judgment without introducing unrealistic extremes.Each factor can vary between very small (near 0) and fairly large (up to \~0.7–0.8), but the average weight remains around 0.20 (1/5). The Gamma(1,1) shape parameter ensures broad variability without producing absurd or hyper-concentrated weights (like one factor getting 99% of total influence).

Each Monte Carlo run represents one hypothetical weighting configuration within those bounds. Monte Carlo error decreases approximately in proportion to the inverse of the square root of the number of simulations:

$$
\text{Error} \propto \frac{1}{\sqrt{N}}
$$

This means that as the number of scenarios (N) increases, the uncertainty in the simulated results shrinks rapidly at first and then levels off:

-   100 runs → approximately 10 % uncertainty\
-   1,000 runs → approximately 3 % uncertainty\
-   5,000 runs → approximately 1.4 % uncertainty

In practice, 5,000 simulations represent a balanced choice—large enough to yield stable percentile and rank estimates, yet small enough to run efficiently on a standard laptop without noticeable slowdown.
:::

### Monte Carlo Variation of Simulated Ranks

@fig-Sensitivity-Analysis summarizes the range of simulated outcomes for each property across all Monte Carlo trials. Each dot represents the property’s median simulated rank, while the horizontal line shows the spread between the 5th and 95th percentile of its rank outcomes. A smaller range means that the property’s position was relatively stable even as the factor weights varied, suggesting its performance is robust under different weighting scenarios. A longer range, on the other hand, indicates that the property’s relative standing fluctuated more widely depending on which factors carried greater influence in a given simulation.

From an appraisal or site-selection standpoint, this visualization conveys the stability of relative rankings under uncertainty. Properties with consistently low (leftward) and narrow intervals can be viewed as strong performers across a broad range of assumptions, implying dependable comparative strength. Conversely, properties with wider or right-shifted intervals are more sensitive to model assumptions, meaning their perceived attractiveness or competitiveness is highly dependent on the specific combination of factor weights.

Taken together, this chart provides a concise summary of ranking volatility, helping analysts gauge not just which sites tend to score best, but how confidently those conclusions hold across multiple possible weighting scenarios.

```{r}
#| label: fig-Sensitivity-Analysis
#| fig-cap: "Sensitivity Analysis: Ranking Stability Under Weight Variation"
#| fig-subcap: "Dot = Median Simulated Rank, Line = 5th–95th Percentile"
#| echo: false
#| message: false
#| warning: false


# ------------------------------------------------------
# 1. Read source file to extract Sale column
# ------------------------------------------------------
src <- read_excel(params$data_file) |> janitor::clean_names()

sale_candidates <- c("sale", "sale_number", "sale_no", "sale_num",
                     "saleid", "sale_id", "sale#", "sale_")

sale_col <- names(src)[tolower(names(src)) %in% sale_candidates]

if (length(sale_col) == 1) {
  sale_lookup <- src |>
    mutate(addr_key = str_trim(tolower(address))) |>
    select(addr_key, Sale = all_of(sale_col))
} else {
  warning("No Sale column found; using Address instead.")
  sale_lookup <- src |>
    mutate(addr_key = str_trim(tolower(address))) |>
    mutate(Sale = address) |>
    select(addr_key, Sale)
}

# ------------------------------------------------------
# 2. Prepare base scores (merge in Sale label)
# ------------------------------------------------------
base_scores <- location_scores %>%
  select(
    Sale = sale,
    Income = income_score,
    Population = density_score,
    Amenities = amenity_score,
    PropertyValue = property_value_score,
    Traffic = traffic_score
  )

# ------------------------------------------------------
# 3. Monte Carlo simulation
# ------------------------------------------------------
set.seed(123)
n_sims <- 500

simulate_rankings <- function() {
  weights <- rgamma(5, shape = 1)
  weights <- weights / sum(weights)
  names(weights) <- c("Income", "Population", "Amenities", "PropertyValue", "Traffic")

  base_scores %>%
    mutate(
      TotalScore = Income * weights["Income"] +
                   Population * weights["Population"] +
                   Amenities * weights["Amenities"] +
                   PropertyValue * weights["PropertyValue"] +
                   Traffic * weights["Traffic"]
    ) %>%
    arrange(desc(TotalScore)) %>%
    mutate(Rank = row_number()) %>%
    select(Sale, Rank)
}

sim_results <- purrr::map_dfr(seq_len(n_sims), ~simulate_rankings(), .id = "Sim")

rank_summary <- sim_results %>%
  group_by(Sale) %>%
  summarise(
    median_rank = median(Rank),
    rank_low = quantile(Rank, 0.05),
    rank_high = quantile(Rank, 0.95),
    .groups = "drop"
  ) %>%
  arrange(median_rank)

# ------------------------------------------------------
# 4. Plot
# ------------------------------------------------------
ggplot(rank_summary, aes(x = reorder(Sale, median_rank))) +
  geom_errorbar(
    aes(ymin = rank_low, ymax = rank_high),
    width = 0.2,
    color = "gray50",
    linewidth = 1
  ) +
  geom_point(
    aes(y = median_rank),
    size = 3,
    color = "steelblue"
  ) +
  coord_flip() +
  scale_y_reverse(breaks = seq(1, nrow(rank_summary), 1)) +
  labs(
    title = "Monte Carlo Variation",
    x = "Sale",
    y = "Simulated Rank (Lower = Better)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )
```
### Simulated Ranking Stability Summary

@tbl-Sensitivity-Analysis-Summary-Table presents the statistical summary of simulated rankings for each sale across all Monte Carlo iterations. The Median Rank column represents the property’s typical position under random variations in factor weights, while the 5th and 95th Percentiles indicate the lower and upper bounds of its likely ranking range. Together, these values convey both the central tendency and the volatility of each property’s simulated performance.

A smaller gap between the 5th and 95th percentiles suggests that the property’s ranking remained relatively stable regardless of how the weights were varied, indicating a higher degree of confidence in its position. Conversely, a wider range implies that the property’s relative standing was more sensitive to the specific weighting assumptions.

In practical terms, this table quantifies the same insights illustrated in the preceding Monte Carlo variation chart. Properties with low median ranks and narrow percentile ranges appear more consistently competitive across numerous simulation scenarios, while those with broader ranges or higher median ranks demonstrate greater uncertainty in their comparative standing. This information helps analysts and appraisers assess both expected performance and model sensitivity, improving the robustness of ranking-based conclusions.

```{r}
#| label: tbl-Sensitivity-Analysis-Summary-Table
#| tbl-cap: "Simulated Ranking Stability Summary"
#| echo: false
#| message: false
#| warning: false


# Format results as a table
rank_summary_table <- rank_summary %>%
  mutate(
    `Median Rank` = round(median_rank, 1),
    `5th Percentile` = round(rank_low, 1),
    `95th Percentile` = round(rank_high, 1)
  ) %>%
  select(Sale, `Median Rank`, `5th Percentile`, `95th Percentile`) %>%
  arrange(`Median Rank`)

# Render formatted table
knitr::kable(
  rank_summary_table,
  digits = 1,
  align = "lccc",
  format = ifelse(knitr::is_latex_output(), "latex", "html")
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    latex_options = "scale_down"
  )

```
### Histogram of Simulated Ranks

@fig-rank-histogram displays the frequency distribution of simulated rankings for each property across all Monte Carlo runs. Each panel corresponds to one sale, and the horizontal axis represents its possible rank outcomes. A tightly clustered distribution suggests that the property’s rank remained consistent despite varying factor weights, indicating a relatively stable competitive position. A broader or multi-peaked distribution reveals greater sensitivity—meaning the property’s relative standing can change noticeably depending on how much weight is assigned to each factor.

From an analytical perspective, this helps identify which properties are most robust to methodological assumptions and which are more dependent on specific weight configurations. In site-selection or appraisal applications, consistent rankings indicate locations that perform well across a wide range of weighting scenarios, while highly variable rankings point to those whose appeal hinges on particular market drivers.

```{r}
#| label: fig-rank-histogram
#| fig-cap: "Histogram of Simulated Ranks by Sale"
#| echo: false
#| message: false
#| warning: false

ggplot(sim_results, aes(x = Rank, fill = Sale)) +
geom_histogram(binwidth = 1, alpha = 0.6, position = "identity") +
facet_wrap(~ Sale, scales = "free_y") +
theme_minimal(base_size = 13) +
labs(
title = "Distribution of Simulated Ranks",
x = "Simulated Rank",
y = "Frequency"
) +
theme(legend.position = "none")
```
### Tornado Diagram of Weight Sensitivity

@fig-tornado presents the distribution of simulated weights assigned to each evaluation factor across all Monte Carlo iterations. Each horizontal boxplot summarizes the spread and central tendency of one factor’s weight. The median line within each box represents the typical relative importance of that factor, while the width of the box reflects how much its weight varied from one simulation to another.

A higher median implies that, on average, the factor contributed more to the overall scores, whereas a lower median suggests it played a smaller role. A narrower spread indicates that its influence was stable across scenarios, while a wider spread denotes greater uncertainty or volatility in its importance. Together, these elements convey how sensitive the overall ranking results are to shifts in the assumed weighting of each variable.

In general, the tornado diagram complements the rank histograms by revealing why certain properties appear more stable or variable—helping analysts trace rank uncertainty back to underlying factors.

```{r}
#| label: fig-tornado
#| fig-cap: "Tornado Diagram — Impact of Weight Variation on Total Score"
#| echo: false

set.seed(123)
n_sims <- 1000
weight_samples <- replicate(n_sims, {
w <- rgamma(5, shape = 1)
w / sum(w)
})

weight_df <- as.data.frame(t(weight_samples))
names(weight_df) <- c("Income", "Population", "Amenities", "PropertyValue", "Traffic")

weight_df_long <- tidyr::pivot_longer(weight_df, everything(), names_to = "Factor", values_to = "Weight")

ggplot(weight_df_long, aes(x = reorder(Factor, Weight, FUN = median), y = Weight)) +
geom_boxplot(fill = "skyblue", color = "gray40") +
coord_flip() +
theme_minimal(base_size = 13) +
labs(
title = "Tornado Diagram of Weight Sensitivity",
x = "Factor",
y = "Simulated Weight Distribution"
)
```

### Radar Chart

The radar chart in @fig-radar illustrates how the top five locations compare across the primary scoring components of the Retail Location Analysis: income, population, amenities, property value, and traffic. Each polygon represents one property’s relative performance, with the outer edge of the chart indicating higher normalized scores. The chart provides a multidimensional view of how well each site performs across the full range of evaluation criteria rather than focusing on a single factor.

From an analytical perspective, a property with a larger and more evenly shaped polygon demonstrates balanced strength across multiple dimensions, while a property whose polygon extends sharply in only one or two directions indicates concentrated advantages but potential weaknesses elsewhere. For example, a site that scores highly on income and traffic but low on amenities or property value may perform well under certain weighting scenarios but less so when broader community attributes are emphasized.

In practice, this visualization helps appraisers and analysts understand why certain sites achieved higher overall scores or ranking stability in the Monte Carlo results. It reveals the underlying trade-offs between economic strength, consumer accessibility, and neighborhood characteristics, allowing decision-makers to assess not just which locations rank best, but how their component profiles differ and where improvement or additional investigation may be warranted.

```{r}
#| label: fig-radar
#| fig-cap: "Component Comparison Radar Chart"
#| echo: false
#| warning: false
#| message: false


# --------------------------------------------------
# 1. Prepare data for radar chart
# --------------------------------------------------

# Select the five scoring components directly from location_scores
radar_data <- location_scores %>%
  select(
    Sale = sale,
    Income = income_score,
    Population = density_score,
    Amenities = amenity_score,
    PropertyValue = property_value_score,
    Traffic = traffic_score,
    Rank = rank
  )

# Keep only top N ranked sites (for clarity)
top_n <- 5
radar_subset <- radar_data %>%
  arrange(Rank) %>%
  slice_head(n = top_n)

# Rescale each variable (except Sale and Rank) to 0–100
radar_scaled <- radar_subset %>%
  mutate(across(c(Income, Population, Amenities, PropertyValue, Traffic), 
                ~ scales::rescale(.x, to = c(0, 100))))

# Convert to long format for radar plotting
radar_long <- radar_scaled %>%
  pivot_longer(
    cols = c(Income, Population, Amenities, PropertyValue, Traffic),
    names_to = "Factor",
    values_to = "Score"
  )

# --------------------------------------------------
# 2. Create radar chart
# --------------------------------------------------

if (knitr::is_html_output()) {

  # Interactive radar (Plotly)
  fig <- plot_ly(type = "scatterpolar", fill = "toself")

  for (site in unique(radar_long$Sale)) {
    df_site <- radar_long %>% filter(Sale == site)
    fig <- fig %>%
      add_trace(
        r = df_site$Score,
        theta = df_site$Factor,
        name = site,
        mode = "lines+markers",
        fill = "toself"
      )
  }

  fig <- fig %>%
    layout(
      polar = list(
        radialaxis = list(visible = TRUE, range = c(0, 100))
      ),
      legend = list(title = list(text = "Sale"), orientation = "v"),
      title = list(
        text = paste0("Component Comparison Radar Chart — Top ", top_n, " Sales"),
        font = list(size = 16)
      )
    )

  fig

} else {

  # Static fallback (ggplot2 for PDF)
  ggplot(radar_long, aes(x = Factor, y = Score, group = Sale, color = Sale)) +
    geom_polygon(fill = NA, linewidth = 0.7, alpha = 0.6) +
    geom_point(size = 2) +
    coord_polar() +
    ylim(0, 100) +
    labs(
      title = paste0("Component Comparison Radar Chart — Top ", top_n, " Sales"),
      y = "Score (0–100)",
      x = NULL
    ) +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.x = element_text(size = 10, face = "bold"),
      plot.title = element_text(face = "bold", hjust = 0.5),
      legend.title = element_text(size = 11, face = "bold"),
      legend.position = "bottom"
    )
}
```
### Ranked Location map

The map in @fig-ranked-location-map provides a spatial summary of the ranked retail locations. Each point represents either the subject property or one of the comparable sales analyzed in the model. The subject property is displayed in red for emphasis, while the comparable sales appear as blue points labeled with their respective ranking positions. This visual layout allows the reader to see not only how each property performed in the scoring analysis but also where those results occur geographically across the market area. Proximity patterns, clustering, or spatial dispersion of higher- or lower-ranked sites can reveal localized strengths and weaknesses that might not be evident from tabular results alone.

```{r}
#| label: fig-ranked-location-map
#| fig-cap: "Ranked Location Map"
#| echo: false
#| message: false
#| warning: false

library(leaflet)
library(dplyr)

# Confirm required data exist
if (!exists("location_scores")) stop("Object 'location_scores' not found — run scoring first.")

# Prepare data
ranked_points <- location_scores %>%
  filter(!is.na(lat) & !is.na(lon)) %>%
  mutate(
    label_text = ifelse(
      sale == "Subject",
      paste0("Subject — Ranked ", rank),
      paste0("Sale ", sale, " — Ranked ", rank)
    )
  )

# Separate subject and comparables
subject_point <- ranked_points %>% filter(sale == "Subject")
sales_points  <- ranked_points %>% filter(sale != "Subject")

# Build leaflet map
leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  # Add comparable sales (simple blue points)
  addCircleMarkers(
    data = sales_points,
    lng = ~lon, lat = ~lat,
    radius = 5,
    stroke = FALSE,
    fillColor = "blue",
    fillOpacity = 1,
    label = ~label_text,
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(12, 0),
      style = list(
        "color" = "black",
        "font-weight" = "bold",
        "font-size" = "12px",
        "text-shadow" = "1px 1px 2px white"
      )
    ),
    group = "Comparable Sales",
    layerId = ~paste0("sale_", sale)
  ) %>%
  
  # Add subject property (red)
  addCircleMarkers(
    data = subject_point,
    lng = ~lon, lat = ~lat,
    radius = 6,
    stroke = FALSE,
    fillColor = "red",
    fillOpacity = 1,
    label = ~label_text,
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(14, 0),
      style = list(
        "color" = "red",
        "font-weight" = "bold",
        "font-size" = "13px",
        "text-shadow" = "1px 1px 2px white"
      )
    ),
    group = "Subject Property",
    layerId = "subject_ranked"
  ) %>%
  
  # Add layer control
  addLayersControl(
    overlayGroups = c("Subject Property", "Comparable Sales"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
  # Add recenter button with automatic recentering
  htmlwidgets::onRender("
    function(el, x) {
      var map = this;
      
      // ============================================
      // CREATE RECENTER BUTTON
      // ============================================
      var RecenterControl = L.Control.extend({
        options: {
          position: 'topleft'
        },
        onAdd: function(map) {
          var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control');
          var button = L.DomUtil.create('a', 'leaflet-control-recenter', container);
          button.innerHTML = '🎯';  // Bullseye emoji;
          button.href = '#';
          button.title = 'Recenter map to visible properties';
          button.style.fontSize = '20px';
          button.style.lineHeight = '30px';
          button.style.width = '30px';
          button.style.height = '30px';
          button.style.textAlign = 'center';
          button.style.textDecoration = 'none';
          button.style.backgroundColor = 'white';
          button.style.cursor = 'pointer';
          
          L.DomEvent.on(button, 'click', function(e) {
            L.DomEvent.preventDefault(e);
            recenterMap();
          });
          
          return container;
        }
      });
      
      map.addControl(new RecenterControl());
      
      // ============================================
      // RECENTER FUNCTION
      // ============================================
      function recenterMap() {
        var bounds = L.latLngBounds();
        var hasVisibleLayers = false;
        
        map.eachLayer(function(layer) {
          if (map.hasLayer(layer) && layer.getBounds) {
            try {
              bounds.extend(layer.getBounds());
              hasVisibleLayers = true;
            } catch(e) {}
          }
          else if (map.hasLayer(layer) && layer.getLatLng) {
            try {
              bounds.extend(layer.getLatLng());
              hasVisibleLayers = true;
            } catch(e) {}
          }
        });
        
        if (hasVisibleLayers && bounds.isValid()) {
          map.fitBounds(bounds, {
            padding: [50, 50],
            maxZoom: 14
          });
        }
      }
      
      // ============================================
      // INITIAL FIT ON LOAD
      // ============================================
      setTimeout(function() {
        recenterMap();
      }, 100);
      
      // ============================================
      // RECENTER WHEN LAYERS TOGGLE
      // ============================================
      map.on('overlayadd overlayremove', function(e) {
        setTimeout(recenterMap, 50);
      });
      
      // ============================================
      // STYLING
      // ============================================
      var style = document.createElement('style');
      style.innerHTML = `
        .leaflet-control-recenter:hover {
          background-color: #f4f4f4;
        }
      `;
      el.appendChild(style);
    }
  ")
```
### Amenity Density Map

The Amenity Density Map in Figure @fig-amenity-density-map expands on the earlier tabular summary by showing the actual spatial distribution of each amenity type within the drive-time boundaries. Each color represents a distinct amenity category, allowing the reader to visualize not only the total number of amenities but also their composition and clustering patterns around each property. Concentrations of similarly colored dots often indicate commercial corridors or service hubs, while mixed clusters suggest diverse neighborhood environments supporting broader retail demand. This map thus provides an intuitive, spatial counterpart to the amenity count analysis, illustrating both density and diversity across trade areas.

```{r}
#| label: build-amenity-points
#| message: false
#| warning: false
#| cache: true

# --------------------------------------------------------------
# Build sf layer of all amenities for each property
# --------------------------------------------------------------

if (!exists("isochrones_sf")) stop("isochrones_sf not found — run isochrone generation first.")
if (!exists("params")) stop("params object missing — ensure your YAML parameters are loaded.")

# Helper function to safely retrieve amenity points
get_amenity_points <- function(sale_id, amenity_type) {
  iso <- isochrones_sf |> filter(sale == sale_id)
  if (nrow(iso) == 0) return(NULL)

  bbox <- st_bbox(iso)
  result <- tryCatch({
    opq(bbox = bbox) |>
      add_osm_feature(key = "amenity", value = amenity_type) |>
      osmdata_sf()
  }, error = function(e) NULL)

  if (is.null(result) || is.null(result$osm_points) || nrow(result$osm_points) == 0) {
    return(NULL)
  }

  pts <- result$osm_points |>
    st_transform(st_crs(iso)) |>
    st_intersection(iso) |>
    mutate(sale = sale_id, amenity_type = amenity_type) |>
    select(osm_id, amenity = amenity_type, sale, geometry)  # standardize columns

  return(pts)
}

# Loop through sales and amenity types
all_amenities_list <- list()
for (s in unique(isochrones_sf$sale)) {
  for (a in params$amenity_types) {
    pts <- get_amenity_points(s, a)
    if (!is.null(pts) && nrow(pts) > 0) {
      all_amenities_list[[paste(s, a, sep = "_")]] <- pts
    }
  }
}

# Harmonize column structure before binding
if (length(all_amenities_list) > 0) {
  amenity_counts_sf <- bind_rows(all_amenities_list)
  amenity_counts_sf <- st_as_sf(amenity_counts_sf)
  amenity_counts_sf <- st_make_valid(amenity_counts_sf)
  message("✅ amenity_counts_sf created with ", nrow(amenity_counts_sf), " points across ",
          length(unique(amenity_counts_sf$sale)), " sales.")
} else {
  warning("No amenity points found — amenity_counts_sf is empty.")
  amenity_counts_sf <- st_sf(osm_id = character(),
                             amenity = character(),
                             sale = character(),
                             geometry = st_sfc(),
                             crs = 4326)
}
```
```{r}
#| label: fig-amenity-density-map
#| fig-cap: "Amenity Density Map"
#| echo: false
#| message: false
#| warning: false

library(leaflet)
library(dplyr)
library(sf)
library(stringr)

# Make sure amenities exist
if (!exists("amenity_counts_sf") || nrow(amenity_counts_sf) == 0) {
  stop("amenity_counts_sf not found — run the amenity retrieval chunk first.")
}

# Ensure consistent CRS
properties_sf <- st_transform(properties_sf, 4326)
isochrones_sf <- st_transform(isochrones_sf, 4326)

# Define base layers
subject_point <- properties_sf %>% filter(sale == "Subject")
sales_points  <- properties_sf %>% filter(sale != "Subject")

# Distinct amenity types
unique_amenities <- sort(unique(amenity_counts_sf$amenity))
pal <- colorFactor(
  palette = RColorBrewer::brewer.pal(length(unique_amenities), "Set2"),
  domain  = unique_amenities
)

# Build neat HTML legend
legend_html <- paste0(
  "<div style='
      background-color: white;
      padding: 8px 12px;
      border-radius: 8px;
      box-shadow: 0 0 5px rgba(0,0,0,0.3);
      font-size: 13px;
      line-height: 1.4;
  '>
   <b>Amenity Type</b><br>",
  paste(
    sprintf(
      "<div style='display:flex; align-items:center; margin-top:2px;'>
         <div style='width:14px; height:14px; background:%s; border:1px solid #555; margin-right:6px;'></div>
         <div>%s</div>
       </div>",
      pal(unique_amenities),
      str_to_title(unique_amenities)
    ),
    collapse = ""
  ),
  "</div>"
)

# Create map
map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  addPolygons(
    data = isochrones_sf,
    fillColor = "#cccccc",
    fillOpacity = 0.1,
    color = "#999999",
    weight = 1,
    group = "Drive-Time Boundaries",
    layerId = ~paste0("iso_", sale)
  ) %>%
  
  addCircleMarkers(
    data = amenity_counts_sf,
    radius = 3,
    stroke = FALSE,
    fillOpacity = 0.8,
    fillColor = ~pal(amenity),
    group = "Amenities",
    label = ~paste0(amenity, " (Sale ", sale, ")"),
    layerId = ~paste0("amenity_", osm_id)
  ) %>%
  
  addCircleMarkers(
    data = sales_points,
    radius = 5,
    stroke = FALSE,
    fillColor = "blue",
    fillOpacity = 1,
    label = ~paste("Sale", sale),
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(10, 0),
      style = list("color" = "black", "font-weight" = "bold", "font-size" = "12px")
    ),
    group = "Comparable Sales",
    layerId = ~paste0("sale_", sale)
  ) %>%
  
  addCircleMarkers(
    data = subject_point,
    radius = 6,
    stroke = FALSE,
    fillColor = "red",
    fillOpacity = 1,
    label = "Subject",
    labelOptions = labelOptions(
      noHide = TRUE,
      direction = "right",
      textOnly = TRUE,
      offset = c(12, 0),
      style = list("color" = "red", "font-weight" = "bold", "font-size" = "13px")
    ),
    group = "Subject Property",
    layerId = "subject_amenity"
  ) %>%
  
  addControl(
    html = legend_html,
    position = "bottomright"
  ) %>%
  
  addLayersControl(
    overlayGroups = c("Amenities", "Drive-Time Boundaries", "Comparable Sales", "Subject Property"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
  # Add recenter button with bullseye icon
  htmlwidgets::onRender("
    function(el, x) {
      var map = this;
      
      // ============================================
      // CREATE RECENTER BUTTON (BULLSEYE ICON)
      // ============================================
      var RecenterControl = L.Control.extend({
        options: {
          position: 'topleft'
        },
        onAdd: function(map) {
          var container = L.DomUtil.create('div', 'leaflet-bar leaflet-control');
          var button = L.DomUtil.create('a', 'leaflet-control-recenter', container);
          button.innerHTML = '🎯';  // Bullseye emoji;
          button.href = '#';
          button.title = 'Recenter map to visible layers';
          button.style.fontSize = '18px';
          button.style.lineHeight = '30px';
          button.style.width = '30px';
          button.style.height = '30px';
          button.style.textAlign = 'center';
          button.style.textDecoration = 'none';
          button.style.backgroundColor = 'white';
          button.style.cursor = 'pointer';
          
          L.DomEvent.on(button, 'click', function(e) {
            L.DomEvent.preventDefault(e);
            recenterMap();
          });
          
          return container;
        }
      });
      
      map.addControl(new RecenterControl());
      
      // ============================================
      // RECENTER FUNCTION
      // ============================================
      function recenterMap() {
        var bounds = L.latLngBounds();
        var hasVisibleLayers = false;
        
        map.eachLayer(function(layer) {
          if (map.hasLayer(layer) && layer.getBounds) {
            try {
              bounds.extend(layer.getBounds());
              hasVisibleLayers = true;
            } catch(e) {}
          }
          else if (map.hasLayer(layer) && layer.getLatLng) {
            try {
              bounds.extend(layer.getLatLng());
              hasVisibleLayers = true;
            } catch(e) {}
          }
        });
        
        if (hasVisibleLayers && bounds.isValid()) {
          map.fitBounds(bounds, {
            padding: [50, 50],
            maxZoom: 14
          });
        }
      }
      
      // ============================================
      // INITIAL FIT ON LOAD
      // ============================================
      setTimeout(function() {
        recenterMap();
      }, 100);
      
      // ============================================
      // RECENTER WHEN LAYERS TOGGLE
      // ============================================
      map.on('overlayadd overlayremove', function(e) {
        setTimeout(recenterMap, 50);
      });
      
      // ============================================
      // STYLING
      // ============================================
      var style = document.createElement('style');
      style.innerHTML = `
        .leaflet-control-recenter:hover {
          background-color: #f4f4f4;
        }
      `;
      el.appendChild(style);
    }
  ")

map

```

# Conclusions

Comparative analysis is defined as "the process by which a value indication is derived in the sales comparison approach. Comparative analysis may employ quantitative or qualitative techniques, either separately or in combination."[^1] The Appraisal of Real Estate outlines examples of the techniques used in quantitative adjustments and qualitative analyses as shown in @tbl-techniques.[^2]

[^1]: *The Dictionary of Real Estate Appraisal*, 7th ed. (Chicago: Appraisal Institute, 2022), s.v. “comparative analysis.”

[^2]: *The Appraisal of Real Estate*, 15th ed. (Chicago: Appraisal Institute, 2020), 362.

```{r}
#| label: tbl-techniques
#| tbl-cap: "Techniques Used in Quantitative and Qualitative Analysis"
#| echo: false
#| warning: false
#| message: false

library(knitr)
library(kableExtra)

# Create the data frame

techniques <- data.frame(
`Quantitative Analysis` = c(
"Paired data analysis (sales and resales of the same or similar properties)",
"Grouped data analysis",
"Secondary data analysis",
"Statistical analysis including graphic analysis and scenario analysis",
"Cost-related adjustments (cost to cure, depreciated cost)",
"Capitalization of income differences",
"Trend analysis"
),
`Qualitative Analysis` = c(
"Relative comparison analysis",
"Ranking analysis",
"Personal interviews",
"",
"",
"",
""
),
stringsAsFactors = FALSE
)

# Render the table cleanly and narrowly

knitr::kable(
techniques,
align = c("l", "l"),
booktabs = TRUE,
col.names = c("Quantitative Analysis", "Qualitative Analysis")
) %>%
kable_styling(
full_width = FALSE,
position = "center",
font_size = 10,
latex_options = "hold_position"
) %>%
column_spec(1, width = "6cm") %>%
column_spec(2, width = "5cm")

```

When adequate data is available, quantitative analysis can be a helpful tool, but often the differences between properties makes it challenging to determine appropriate adjustments with a level of certainty. In these cases qualitative analysis or ranking can be helpful."[^3] Qualitative analysis recognizes the inefficiencies of real estate markets and the difficulty of expressing adjustments with mathematical precision."

[^3]: *Appraisal of Real Estate*, 15th ed. , 376.

## Executive Summary Table

@tbl-Executive-Summary consolidates the results of the Retail Location Analysis into a concise summary of comparative performance. Each property (including the subject) is ranked according to its Total Score, which integrates multiple weighted factors such as income, population, amenities, property value, and traffic conditions within the specified drive-time boundary.

The Rank column reflects the overall ordering from best to least favorable under the baseline weighting scenario, while the Total Score expresses each location’s composite index value on a standardized scale. Drive Time denotes the consistent spatial boundary used for comparison, and Price per Square Foot provides a market context for relative value.

Taken together, this table provides a high-level snapshot of how each location performs under the modeled conditions, allowing for quick identification of competitive strengths and potential trade-offs between market positioning and pricing.

```{r}
#| label: tbl-Executive-Summary
#| tbl-cap: "Executive Summary — Retail Location Comparison (Drive-Time Isochrones)"
#| echo: false
#| warning: false
#| message: false

# Read Excel source file to retrieve Sale and Price/SF
src <- read_excel(params$data_file) |> clean_names()

sale_candidates <- c("sale", "sale_number", "sale_no", "sale_num",
                     "saleid", "sale_id", "sale#", "sale_")
sale_col <- names(src)[tolower(names(src)) %in% sale_candidates]

price_candidates <- c("price_sf", "price_per_sf", "price_psf", "price_sqft", "price_sq_ft")
price_col <- names(src)[tolower(names(src)) %in% price_candidates]

lat_col <- intersect(names(src), c("lat", "latitude"))[1]
lon_col <- intersect(names(src), c("lon", "longitude", "long", "x"))[1]

# Build lookup table from Excel
if (length(sale_col) == 1 && length(price_col) == 1 && !is.na(lat_col) && !is.na(lon_col)) {
  src_lookup <- src |>
    mutate(
      lat_join = coalesce(!!sym(lat_col), NA_real_),
      lon_join = coalesce(!!sym(lon_col), NA_real_)
    ) |>
    select(Sale = all_of(sale_col[1]),
           Price_SF = all_of(price_col[1]),
           lat_join, lon_join)
} else {
  stop("Could not locate Sale, Price/SF, or coordinate columns in Excel file.")
}

# Join Excel Sale + Price/SF data back into location_scores
summary_data <- location_scores |>
  mutate(
    lat_join = round(as.numeric(lat), 4),
    lon_join = round(as.numeric(lon), 4)
  ) |>
  left_join(
    src_lookup |>
      mutate(
        lat_join = round(lat_join, 4),
        lon_join = round(lon_join, 4)
      ),
    by = c("lat_join", "lon_join")
  ) |>
  mutate(
    Rank = rank,
    `Total Score` = round(total_score, 1),
    `Drive Time` = paste0(selected_drive_time, " min"),
    `Price/SF` = case_when(
      is.na(Price_SF) ~ NA_character_,
      str_detect(as.character(Price_SF), "^\\$") ~ as.character(Price_SF),
      TRUE ~ dollar(as.numeric(Price_SF), accuracy = 1)
    )
  ) |>
  arrange(Rank) |>
  select(Rank, Sale, `Total Score`, `Drive Time`, `Price/SF`)

# Render table
tbl <- kable(
  summary_data,
  digits = 1,
  format = ifelse(knitr::is_latex_output(), "latex", "html"),
  booktabs = TRUE,
  longtable = FALSE,
  align = c("c", "c", "c", "c", "c")
) |>
  kable_styling(
    full_width = FALSE,
    position = "center",
    latex_options = c("hold_position", "scale_down")
  ) |>
  row_spec(0, bold = TRUE)

if (knitr::is_latex_output()) {
  knitr::asis_output(tbl)
} else {
  tbl
}
```



## Session Information

```{r}
sessionInfo()
```
