Merge pull request #43 from cal-itp/stm-report

STM report cont

Author: natam1

Makefile

@@ -0,0 +1,20 @@
+
+add_precommit:
+	pip install pre-commit
+	pre-commit install
+
+# install_env:
+#	pip install uv && uv sync --all-groups
+#	make add_precommit
+
+install_env:
+	pip install uv && uv sync --all-groups
+	make add_precommit
+
+
+#uv_setup_project:
+#	pip install uv
+#	uv init
+#	uv add shared-utils
+#   uv add calitp-data-analysis
+#	uv lock

_shared_utils/pyproject.toml

@@ -0,0 +1,54 @@
+[project]
+name = "shared-utils"
+version = "5.0.5"
+requires-python = ">=3.11.0, <3.12.0"
+dependencies = [
+    "calitp-data-analysis",
+    "altair==6.0.0",
+    "altair-transform==0.2.0",
+    "branca>=0.7.0",
+    "dask>=2024.10.0,<2025.4.0",
+    "dask-geopandas>=0.5.0",
+    "folium>=0.20.0,<0.21.0",
+    "gcsfs==2024.5.0",
+    "great_tables==0.20.0",
+    "intake==0.6.4",
+    "intake-geopandas>=0.3.0",
+    "intake-parquet",
+    "itables==2.7.3",
+    "loguru>=0.6.0,<0.7",
+    "mapclassify>=2.8.0,<3",
+    "matplotlib>=3.10.5,<4.0.0",
+    "movingpandas==0.22.4",
+    "msgpack>=1.1.2",
+    "numba>=0.62.1,<0.63.0",
+    "numpy>=1.26.4,<2.0.0",
+    "omegaconf==2.3.0",
+    "openpyxl>=3.1.0",
+    "polars==1.22.0",
+    "python-slugify>=8.0.0",
+    "pyyaml>=6.0",
+    "quarto==0.1.0",
+    "quarto-cli==1.6.40; sys_platform != 'win32'",
+    "scipy>=1.11.0",
+    "seaborn>=0.11.2,<0.12.0",
+    "stonesoup==1.8",
+    "tqdm>=4.66.0",
+    "xmltodict>=0.13.0,<0.14",
+    "pyaml>=21.10.1",
+]
+
+# [dependency-groups]
+# test = [
+#     "pytest>=9.0.2,<10.0.0",
+#     "pytest-mock>=3.15.1,<4.0.0",
+#     "pytest-recording>=0.13.4,<0.14.0",
+#     "pytest-unordered>=0.7.0,<0.8.0",
+# ]
+
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[tool.hatch.build.targets.wheel]
+packages = ["shared_utils"]

_shared_utils/shared_utils/arcgis_query.py

@@ -0,0 +1,137 @@
+"""
+Import ESRI feature layers from open data portal sources.
+Export into shared GCS folder to use across projects.
+
+https://github.com/cagov/caldata-mdsa-caltrans-pems/blob/main/jobs/utils/geo.py
+"""
+
+from urllib.parse import parse_qsl, urlencode, urlparse, urlunparse
+
+import geopandas as gpd
+import pandas as pd
+
+from calitp_data_analysis import utils
+from calitp_data_analysis.sql import to_snakecase
+
+SHARED_GCS = "gs://calitp-analytics-data/data-analyses/shared_data/"
+
+COUNTY_POLYGONS_URL = "https://services1.arcgis.com/jUJYIo9tSA7EHvfZ/arcgis/rest/services/California_County_Boundaries/FeatureServer/0/query"
+
+CALTRANS_DISTRICTS_URL = (
+    "https://caltrans-gis.dot.ca.gov/arcgis/rest/services/CHboundary/District_Tiger_Lines/FeatureServer/0/query/"
+)
+
+LEGISLATIVE_BASE = "https://services3.arcgis.com/fdvHcZVgB2QSRNkL/arcgis/rest/services/Legislative/FeatureServer/"
+
+LEGISLATIVE_DICT = {
+    "ca_assembly_districts": f"{LEGISLATIVE_BASE}0/query/",
+    "ca_senate_districts": f"{LEGISLATIVE_BASE}1/query/",
+    "ca_congressional_districts": f"{LEGISLATIVE_BASE}2/query/",
+}
+
+CALTRANS_BASE = "https://caltrans-gis.dot.ca.gov/arcgis/rest/services/CHhighway/"
+CRS_FUNCTIONAL_CLASSICIATION_URL = f"{CALTRANS_BASE}CRS_Functional_Classification/FeatureServer/0/query/"
+SHN_LINES_URL = f"{CALTRANS_BASE}SHN_Lines/FeatureServer/0/query/"
+SHN_POSTMILES_URL = f"{CALTRANS_BASE}SHN_Postmiles_Tenth/FeatureServer/0/query/"
+
+
+def gdf_from_esri_feature_service(url):
+    """
+    Load an Esri Feature Service to a GeoDataFrame.
+
+    Given a URL to an Esri Feature Service, download the features
+    as GeoJSON, and put them into a GeoDataFrame.
+    """
+    parsed = urlparse(url)
+
+    # Ensure we are using the query endpoint of the feature service
+    if not parsed.path.endswith("/query"):
+        parsed = parsed._replace(path=parsed.path + "/query")
+
+    # Keep grabbing data using the resultOffset until there is no more left
+    offset = 0
+    gdfs = []
+    while True:
+        queries = dict(parse_qsl(parsed.query))
+        queries.update(
+            {
+                "where": "1=1",  # Ensure all rows
+                "f": "geojson",  # Ensure GeoJSON
+                "outFields": "*",  # Ensure all columns
+                "resultOffset": str(offset),  # offset the start
+                "returnGeometry": "true",  # Yes we want geometries
+            }
+        )
+        offset_url = urlunparse(parsed._replace(query=urlencode(queries)))
+
+        gdf = gpd.read_file(offset_url, driver="GeoJSON")
+        if len(gdf) == 0:
+            break
+
+        gdfs.append(gdf)
+        offset += len(gdf)
+    return pd.concat(gdfs).reset_index(drop=True)
+
+
+def combine_legislative_districts(assembly_districts_url: str, senate_districts_url: str) -> gpd.GeoDataFrame:
+    """
+    Create a combined assembly district and senate districts
+    gdf.
+    """
+    assembly_districts = gdf_from_esri_feature_service(assembly_districts_url)
+    senate_districts = gdf_from_esri_feature_service(senate_districts_url)
+
+    gdf = pd.concat(
+        [
+            assembly_districts[["AssemblyDistrictLabel", "geometry"]].rename(
+                columns={"AssemblyDistrictLabel": "legislative_district"}
+            ),
+            senate_districts[["SenateDistrictLabel", "geometry"]].rename(
+                columns={"SenateDistrictLabel": "legislative_district"}
+            ),
+        ],
+        axis=0,
+        ignore_index=True,
+    )
+
+    return gdf
+
+
+def exclude_columns(
+    gdf: gpd.GeoDataFrame, list_of_cols: list = ["objectid", "shape__area", "shape__length"]
+) -> gpd.GeoDataFrame:
+    """
+    Drop a couple of columns that tend to show up for ESRI.
+    """
+    for c in list_of_cols:
+        if c in gdf.columns:
+            gdf = gdf.drop(columns=c)
+
+    return gdf
+
+
+if __name__ == "__main__":
+    esri_datasets = {
+        "ca_county": COUNTY_POLYGONS_URL,
+        "caltrans_districts": CALTRANS_DISTRICTS_URL,
+        "ca_congressional_districts": LEGISLATIVE_DICT["ca_congressional_districts"],
+        "state_highway_network_raw": SHN_LINES_URL,
+        "state_highway_network_postmiles": SHN_POSTMILES_URL,
+        "public_road_functional_classification": CRS_FUNCTIONAL_CLASSICIATION_URL,
+    }
+    for dataset_name, url in esri_datasets.items():
+        print(dataset_name)
+        gdf = gdf_from_esri_feature_service(url)
+        gdf = gdf.pipe(to_snakecase).pipe(exclude_columns)
+
+        print(gdf.crs)
+        print(gdf.shape)
+        utils.geoparquet_gcs_export(gdf, SHARED_GCS, dataset_name)
+        del gdf
+
+    legislative_districts_gdf = combine_legislative_districts(
+        LEGISLATIVE_DICT["ca_assembly_districts"],
+        LEGISLATIVE_DICT["ca_senate_districts"],
+    )
+
+    utils.geoparquet_gcs_export(legislative_districts_gdf, SHARED_GCS, "legislative_districts")

_shared_utils/shared_utils/bq_utils.py

@@ -0,0 +1,67 @@
+"""
+pandas_gbq utils to download Big Query tables
+"""
+
+from typing import Literal
+
+import geopandas as gpd
+import google.auth
+import pandas as pd
+import pandas_gbq
+from shared_utils import geo_utils
+
+credentials, project = google.auth.default()
+
+
+def basic_sql_query(project_name: str, dataset_name: str, table_name: str) -> str:
+    """
+    Set up the basic sql query needed, which is the entire table.
+    """
+    sql_query = f"SELECT * FROM  `{project_name}`.`{dataset_name}`.`{table_name}`"
+
+    return sql_query
+
+
+def add_sql_date_filter(date_col: str, start_date: str, end_date: str) -> str:
+    """
+    Add a where condition to filter by date, coerce the dates so sql_query is read correctly.
+    """
+    where_condition = f"WHERE {date_col} >= DATE('{start_date}') AND {date_col} <= DATE('{end_date}')"
+
+    return where_condition
+
+
+def download_table(
+    project_name: str = "cal-itp-data-infra",
+    dataset_name: str = "mart_gtfs",
+    table_name: str = "",
+    date_col: Literal["service_date", "month_first_day", None] = "",
+    start_date: str = "",
+    end_date: str = "",
+    geom_col: str = None,
+    geom_type: Literal["point", "line"] = None,
+) -> Literal[pd.DataFrame, gpd.GeoDataFrame]:
+    """
+    Set up a basic query and use pandas_gbq to import.
+    Coerce datetime column and convert to gdf if needed.
+    """
+    basic_query = basic_sql_query(project_name, dataset_name, table_name)
+    where_condition = add_sql_date_filter(date_col, start_date, end_date)
+    sql_query_statement = f"{basic_query} {where_condition}"
+
+    if date_col is None:
+        df = pandas_gbq.read_gbq(basic_query, project_id=project_name, dialect="standard", credentials=credentials)
+
+        print(f"query: {basic_query}")
+
+    if date_col is not None:
+        df = pandas_gbq.read_gbq(
+            sql_query_statement, project_id=project_name, dialect="standard", credentials=credentials
+        ).astype({date_col: "datetime64[ns]"})
+
+        print(f"query: {sql_query_statement}")
+
+    if geom_col is not None:
+        df = geo_utils.convert_to_gdf(df, geom_col, geom_type)
+
+    return df

_shared_utils/shared_utils/catalog_utils.py

@@ -0,0 +1,45 @@
+"""
+Functions for opening yaml catalogs in shared_utils.
+"""
+
+from pathlib import Path
+from typing import Literal
+
+import intake
+from omegaconf import OmegaConf  # this is yaml parser
+
+shared_utils_directory = "data-analyses/_shared_utils/shared_utils/"
+
+
+def get_catalog_file(catalog_name, home_path=Path.home(), current_path=Path.cwd()):
+    filename = f"{shared_utils_directory}{catalog_name}.yml"
+    parent_directory = current_path
+
+    if home_path not in current_path.parents:
+        raise RuntimeError("The data-analyses repo should be located in your home directory.")
+
+    while True:
+        test_path = parent_directory.joinpath(filename)
+
+        if test_path.is_file():
+            return test_path
+
+        if parent_directory == home_path:
+            raise FileNotFoundError(f"No such catalog file found: {filename}")
+
+        parent_directory = parent_directory.parent
+
+
+def get_catalog(catalog_name: Literal["shared_data_catalog", "gtfs_analytics_data"]) -> Path:
+    """
+    Grab either the shared_data_catalog (uses intake driver) or
+    gtfs_analytics_data catalog (uses OmegaConf yaml parser).
+
+    """
+    catalog_path = get_catalog_file(catalog_name)
+
+    if catalog_name == "gtfs_analytics_data":
+        return OmegaConf.load(catalog_path)
+
+    else:
+        return intake.open_catalog(catalog_path)