<jupyter_start><jupyter_text>One of the things we've struggled with for a long time is having a persistent online reference point for the different EPA Ecoregion datasets that we use in biogeography. Data files are online and released via [EPA's web site](https://www.epa.gov/eco-research/ecoregions), but that's been problematic at times when EPA has either taken things down for a while or moved them. Ideally, we would have some kind of persistent ID that we could always reference, that ID would be resolvable with code so we could reference it and get additional information for software systems, and we would have a place to rally value-added properties and linkages over time. This seemed like another reasonable use case for WikiData in some ongoing research I'm doing on uses of that system for other purposes.In looking through WikiData any hint of these records already being online, I found items for [Level III Ecoregions](https://www.wikidata.org/wiki/Q52111338) and [Level IV Ecoregions](https://www.wikidata.org/wiki/Q52111409). However, there are only 20 of 182 level 3 ecoregions actually created as instances of that item and 44 level 4 ecoregions. They were created by the same user some time ago. I reached out for clarification on intent to [User:TimK MSI](https://www.wikidata.org/wiki/User:TimK_MSI) who came up with a fairly reasonable start to a data model with items like [Rocky Mountains](https://www.wikidata.org/wiki/Q67197232). This includes some variation but mostly has the following elements:* Link to the country (all US as the focus has been on the CONUS mapping)* Links to US States the ecoregion intersects* Representative point location as coordinates* Some have "has part" links to Level 4 ecoregion items* Some have relationships to similarly named items showing a connection to things like physiographic featuresOne issue I see is that the current items use a series ordinal qualifier on the Level III Ecoregion valued class property as a way of capturing an ID that corresponds to the numbering system specific to the CONUS derivation of the North American Ecoregions dataset. A better approach would be to formalize the actual identifier properties and assign the IDs there. I will work through that process as it will be fun to see how that plays out in the WikiData community, but that will be a bit more involved. Having the identifiers on board the WikiData items themselves will allow existing information systems using things like "US_L3CODE" as a key to readily query WikiData with known information.Data prep in this case is pretty straightforward but will involve a number of steps to translate all of the relevant data into appropriate WikiData items and statements. To be fully complete, we should incorporate all ecoregion concepts as defined by Omernik and others in the original work from Level I through Level IV and harmonizing the CONUS derivative of L3 ecoregions with the North American ecoregions by linking their IDs. For this purpose, we are going to need to get all of the individual source data for level 1-4 ecoregions and stitch everything together.I built out a separate process to develop a political boundaries dataset for all three countries covering the North America ecoregions. This is loaded to a PostGIS instance for spatial processing to obtain the intersections between ecoregion boundaries and political boundaries of interest.I looked into two other things I would have liked to do with these data: getting the full multipolygon geometry online somewhere for each ecoregion boundary and loading the ecoregion description information online that is currently locked up in a downloadable PDF. Unfortunately, both of those objectives are a little problematic.Sticking within the Wikiverse and according to documentation, we should be able to load the geometry to Wikimedia Commons items and reference them with the geoshape property. However, after looking into where and how this has been done in other items, I really didn't see any added value to that approach. Any programming against the system would have to go through some pretty circuitous routes to put everything together for use in analysis or visualization. What we would want is some kind of immediate access to a geospatial service from the WikiData item that would return the geospatial part of the data in a variety of formats. Unfortunately, I wasn't able to track down anything close to an authoritative, stable source for that, even from our own geospatial servers in USGS. An interesting phenomenon seen online are numerous cases where people have clipped out an ecoregion subset of one kind or another at State boundaries, probably for ecosystem context and analysis within the State. I've yet to find one of those that shows a careful provenance trace to where the information comes from or how it was processed.I should be able to do something interesting with the ecoregion descriptions, but I need to think about it a while longer. The original source is in a PDF provided from the same source as the data and released in 2011 by the are a little bit in a government report type of release from the Commission for Environmental Cooperation. That document has a Creative Commons-No Derivatives license, which essentially prohibits taking pieces of it out and putting them somewhere else. However, EPA also released a Word version of the basic descriptive content that should be completely public domain, and I could conceivably build some code to scrape the individual descriptions and create Wikimedia Commons articles for them as a reference point. That would give us a persistent identity and the content to work with. But there are a number of Wikipedia articles that mostly describe some of the same basic ecosystems, albeit from a different perspective with Omernik's work just one source. I haven't figured out the most elegant way to link these various artifacts together just yet.<jupyter_code>%matplotlib inline
from functions import get_source_df, lookup_wd
import os
import copy
from wikidataintegrator import wdi_helpers, wdi_core
import pickle

import geopandas as gpd
import pandas as pd
from sqlalchemy import create_engine
from geoalchemy2 import Geometry, WKTElement<jupyter_output><empty_output><jupyter_text>I built all workflow steps here into a suite of functions that can be composed in a variety of ways but generally need to operate in sequence.<jupyter_code>def get_source_data():
    na_er1 = get_source_df(
        source_url="ftp://newftp.epa.gov/EPADataCommons/ORD/Ecoregions/cec_na/na_cec_eco_l1.zip",
        add_columns=[('source', 'NA_CEC_Eco_Level1')],
        target_file_name="NA_CEC_Eco_Level1.shp"
    )

    na_er2 = get_source_df(
        source_url="ftp://newftp.epa.gov/EPADataCommons/ORD/Ecoregions/cec_na/na_cec_eco_l2.zip",
        add_columns=[('source', 'NA_CEC_Eco_Level2')],
        target_file_name="NA_CEC_Eco_Level2.shp"
    )

    na_er3 = get_source_df(
        source_url="ftp://newftp.epa.gov/EPADataCommons/ORD/Ecoregions/cec_na/NA_CEC_Eco_Level3.zip",
        add_columns=[('source', 'NA_CEC_Eco_Level3')]
    )

    conus_er3 = get_source_df(
        source_url="ftp://newftp.epa.gov/EPADataCommons/ORD/Ecoregions/us/us_eco_l3.zip",
        add_columns=[('source', 'US_Eco_Level3')]
    )

    conus_er4 = get_source_df(
        source_url="ftp://newftp.epa.gov/EPADataCommons/ORD/Ecoregions/us/us_eco_l4.zip",
        add_columns=[('source', 'US_Eco_Level4')],
        target_file_name="us_eco_l4_no_st.shp"
    )
    
    na_er1 = na_er1.to_crs({'init': 'epsg:4326'})
    na_er2 = na_er2.to_crs({'init': 'epsg:4326'})
    na_er3 = na_er3.to_crs({'init': 'epsg:4326'})
    conus_er3 = conus_er3.to_crs({'init': 'epsg:4326'})
    conus_er4 = conus_er4.to_crs({'init': 'epsg:4326'})
    
    return na_er1, na_er2, na_er3, conus_er3, conus_er4

def clean_synth_df(na_er1, na_er2, na_er3, conus_er3, conus_er4):
    na_er1["contextual_identifier"] = na_er1.apply(lambda row: f"NA_L1CODE:{row.NA_L1CODE}", axis=1)
    na_er1["common_name"] = na_er1.apply(lambda row: row.NA_L1NAME, axis=1)
    na_er1.drop(
        [c for c in list(na_er1.columns) if c not in ["source","contextual_identifier","common_name","geometry"]], 
        axis=1, 
        inplace=True
    )

    na_er2["contextual_identifier"] = na_er2.apply(lambda row: f"NA_L2CODE:{row.NA_L2CODE}", axis=1)
    na_er2["common_name"] = na_er2.apply(lambda row: row.NA_L2NAME, axis=1)
    na_er2["part_of_identifiers"] = na_er2.apply(lambda row: f"NA_L1CODE:{row.NA_L1CODE}", axis=1)
    na_er2.drop(
        [c for c in list(na_er2.columns) if c not in ["source","contextual_identifier","part_of_identifiers","common_name","geometry"]], 
        axis=1, 
        inplace=True
    )

    na_er3["contextual_identifier"] = na_er3.apply(lambda row: f"NA_L3CODE:{row.NA_L3CODE}", axis=1)
    na_er3["common_name"] = na_er3.apply(lambda row: row.NA_L3NAME, axis=1)
    na_er3["part_of_identifiers"] = na_er3.apply(lambda row: f"NA_L1CODE:{row.NA_L1CODE};NA_L2CODE:{row.NA_L2CODE}", axis=1)
    na_er3.drop(
        [c for c in list(na_er3.columns) if c not in ["source","contextual_identifier","part_of_identifiers","common_name","geometry"]], 
        axis=1, 
        inplace=True
    )

    conus_er3["contextual_identifier"] = conus_er3.apply(lambda row: f"US_L3CODE:{row.US_L3CODE}", axis=1)
    conus_er3["common_name"] = conus_er3.apply(lambda row: row.US_L3NAME, axis=1)
    conus_er3["part_of_identifiers"] = conus_er3.apply(lambda row: f"NA_L1CODE:{row.NA_L1CODE};NA_L2CODE:{row.NA_L2CODE};NA_L3CODE:{row.NA_L3CODE}", axis=1)
    conus_er3.drop(
        [c for c in list(conus_er3.columns) if c not in ["source","contextual_identifier","part_of_identifiers","common_name","geometry"]], 
        axis=1, 
        inplace=True
    )

    conus_er4["contextual_identifier"] = conus_er4.apply(lambda row: f"US_L4CODE:{row.US_L4CODE}", axis=1)
    conus_er4["common_name"] = conus_er4.apply(lambda row: row.US_L4NAME, axis=1)
    conus_er4["part_of_identifiers"] = conus_er4.apply(lambda row: f"NA_L1CODE:{row.NA_L1CODE};NA_L2CODE:{row.NA_L2CODE};NA_L3CODE:{row.NA_L3CODE};US_L3CODE:{row.US_L3CODE}", axis=1)
    conus_er4.drop(
        [c for c in list(conus_er4.columns) if c not in ["source","contextual_identifier","part_of_identifiers","common_name","geometry"]], 
        axis=1, 
        inplace=True
    )

    return na_er1, na_er2, na_er3, conus_er3, conus_er4

def combine_source_data(na_er1, na_er2, na_er3, conus_er3, conus_er4):
    ecoregions = gpd.GeoDataFrame(
        pd.concat(
            [
                na_er1, 
                na_er2, 
                na_er3, 
                conus_er3,
                conus_er4
            ], 
            ignore_index=True, 
            sort=False
        )
    )
    ecoregions.crs = na_er1.crs
    
    return ecoregions

def combine_geometries(combined_source):
    return combined_source.dissolve(by='contextual_identifier')

def prepare_wikidata_info(combined_source):
    combined_source.reset_index(level=0, inplace=True)

    combined_source["x"] = combined_source.centroid.x
    combined_source["y"] = combined_source.centroid.y

    combined_source['common_name'] = combined_source['common_name'].str.title()
    combined_source["wikidata_id"] = combined_source.loc[
        combined_source.source.isin(["US_Eco_Level3","US_Eco_Level4"])
    ].apply(lambda row: lookup_wd(row["contextual_identifier"]), axis=1)
    combined_source = combined_source.where(pd.notnull(combined_source), None)
    
    return combined_source
    
def cache_source(combined_source):
    combined_source.to_file("data_cache/ecoregions_integrated.gpkg", driver="GPKG")
    
def load_cached_source():
    return gpd.read_file("data_cache/ecoregions_integrated.gpkg")

def prep_source_for_pg(combined_dataset):
    dataset_multi = combined_dataset.apply(copy.deepcopy)
    dataset_poly = combined_dataset.apply(copy.deepcopy)
    dataset_multi.drop(
        dataset_multi[dataset_multi.geometry.geom_type == "Polygon"].index, 
        inplace=True
    )
    dataset_poly.drop(
        dataset_poly[dataset_poly.geometry.geom_type == "MultiPolygon"].index, 
        inplace=True
    )

    dataset_multi['geom'] = dataset_multi['geometry'].apply(lambda x: WKTElement(x.wkt, srid=4326))
    dataset_poly['geom'] = dataset_poly['geometry'].apply(lambda x: WKTElement(x.wkt, srid=4326))

    dataset_multi.drop('geometry', 1, inplace=True)
    dataset_poly.drop('geometry', 1, inplace=True)
    
    return dataset_multi, dataset_poly

def get_db_engine():
    db_connection_url = f"postgres://{os.environ.get('PG_USER')}:{os.environ.get('PG_PASS')}@{os.environ.get('PG_HOST')}:{os.environ.get('PG_PORT')}/ecoregions";
    return create_engine(db_connection_url)

def load_to_pg(dataset_multi, dataset_poly):
    db_engine = get_db_engine()
    
    dataset_poly.to_sql(
        "ecoregions_p", 
        db_engine, 
        if_exists='replace', 
        index=True,
        dtype={'geom': Geometry('POLYGON', srid=4326)}
    )

    dataset_multi.to_sql(
        "ecoregions_m", 
        db_engine, 
        if_exists='replace', 
        index=True,
        dtype={'geom': Geometry('MULTIPOLYGON', srid=4326)}
    )
    
def get_intersections():
    return pd.read_sql("SELECT * FROM intersections", get_db_engine())

def boundary_intersects_list(contextual_identifier, intersect_type="admin_id"):
    intersecting_ids = df_intersections.loc[df_intersections.contextual_identifier == contextual_identifier][intersect_type].to_list()
    intersecting_ids = list(set(intersecting_ids))
    return intersecting_ids

def recombine_ecoregions():
    db_engine = get_db_engine()
    
    q_ecoregions = '''SELECT 
        contextual_identifier,
        wikidata_id,
        source, 
        common_name, 
        part_of_identifiers, 
        x, 
        y 
        FROM ecoregions_m
        UNION
        SELECT 
        contextual_identifier, 
        wikidata_id,
        source, 
        common_name, 
        part_of_identifiers, 
        x, 
        y 
        FROM ecoregions_p
        '''

    df_ecoregions = pd.read_sql(q_ecoregions, db_engine)
    df_ecoregions['part_of_identifiers'] = df_ecoregions['part_of_identifiers'].apply(lambda x: x.split(';') if x is not None else x)
    df_ecoregions['admin_intersects'] = df_ecoregions['contextual_identifier'].apply(lambda x: boundary_intersects_list(x, "admin_id"))
    df_ecoregions['country_intersects'] = df_ecoregions['contextual_identifier'].apply(lambda x: boundary_intersects_list(x, "country_id"))
    
    return df_ecoregions

def cache_final_output(final_data):
    outfile = open("data_cache/ecoregions.pkl", "wb")
    pickle.dump(final_data.to_dict(orient="records"), outfile)
    outfile.close()<jupyter_output><empty_output><jupyter_text>I built some helper functions to deal with common chores in this workflow, including functions to retrieve source data from FTP or HTTP URLs, run some common transformation steps, and return geodataframes for processing. The following function retrieves all 5 ecoregion source datasets as shapefiles and runs basic transformation steps.<jupyter_code>%%time
na_er1, na_er2, na_er3, conus_er3, conus_er4 = get_source_data()<jupyter_output>CPU times: user 26 s, sys: 1.23 s, total: 27.2 s
Wall time: 28.2 s<jupyter_text>In this function, I simplify and standardize the data schemas of each ecoregion source dataset so to make things simpler for further processing into WikiData items. One of the main things here is setting up each level of the inherently hierarchical system with the ecoregion identifiers that will eventually become links to the related items in WikiData.<jupyter_code>%%time
na_er1, na_er2, na_er3, conus_er3, conus_er4 = clean_synth_df(na_er1, na_er2, na_er3, conus_er3, conus_er4)<jupyter_output>CPU times: user 546 ms, sys: 7.78 ms, total: 554 ms
Wall time: 563 ms<jupyter_text>In this codeblock, I concatenate all ecoregion source data together into one dataframe now that everything has the same schema.<jupyter_code>%%time
ecoregions = combine_source_data(na_er1, na_er2, na_er3, conus_er3, conus_er4)<jupyter_output>CPU times: user 7.68 ms, sys: 1.5 ms, total: 9.18 ms
Wall time: 8.95 ms<jupyter_text>In subsequent steps, I need to get a representational point for each ecoregion unit and run intersections with administrative boundaries. To make this a bit simpler in providing values grouped/centered on each distinct ecoregion unit, I run a dissolve and group by the contextual_identifier created earlier.<jupyter_code>%%time
ecoregions = combine_geometries(ecoregions)<jupyter_output>CPU times: user 2min 44s, sys: 6.12 s, total: 2min 50s
Wall time: 2min 58s<jupyter_text>In this step, I introduce x and y coordinates for a representational point to add to the WikiData items, proper case the ecoregion name for the WikiData label, and run a lookup on a previously built WikiData reference set where I pulled in all current instances of Level 3 and 4 CONUS ecoregions.<jupyter_code>%%time
ecoregions = prepare_wikidata_info(ecoregions)<jupyter_output>CPU times: user 5.4 s, sys: 236 ms, total: 5.64 s
Wall time: 6.2 s<jupyter_text>At this stage, we have a fairly reasonable source dataset that has everything in it except for the computationally intensive part of intersecting with administrative boundaries. To avoid having to run all of the above processes again, I cache out an intermediary product that can be reloaded as needed.<jupyter_code>ecoregions.head()
%%time
cache_source(ecoregions)<jupyter_output>CPU times: user 33.7 s, sys: 4.27 s, total: 38 s
Wall time: 41.1 s<jupyter_text>In this step, I load the two datasets of like geometry to PostGIS tables for spatial processing to find intersections with administrative boundaries.<jupyter_code>%%time
dataset_multi, dataset_poly = prep_source_for_pg(ecoregions)
%%time
load_to_pg(dataset_multi, dataset_poly)<jupyter_output>CPU times: user 2.57 s, sys: 1.51 s, total: 4.08 s
Wall time: 1min 7s<jupyter_text>To finish this part of things out for further processing into WikiData items, I used processes directly in PostGIS to have it handle the more intensive spatial processing. This involved creating a materialized view containing the names and identifiers of intersected boundaries that can be used to generate WikiData links (statements/claims) to relevant items. This was handled with the following SQL.    SELECT b.name AS boundary_name,     b.wikidata_id AS wikidata_id,    e.contextual_identifier    FROM boundaries_m AS b, ecoregions_m AS e    WHERE ST_Intersects(b.geom, e.geom)    UNION    SELECT b.name AS boundary_name,     b.wikidata_id AS wikidata_id,    e.contextual_identifier    FROM boundaries_m AS b, ecoregions_p AS e    WHERE ST_Intersects(b.geom, e.geom)    UNION    SELECT b.name AS boundary_name,     b.wikidata_id AS wikidata_id,    e.contextual_identifier    FROM boundaries_p AS b, ecoregions_m AS e    WHERE ST_Intersects(b.geom, e.geom)    UNION    SELECT b.name AS boundary_name,     b.wikidata_id AS wikidata_id,    e.contextual_identifier    FROM boundaries_p AS b, ecoregions_p AS e    WHERE ST_Intersects(b.geom, e.geom)Note: I do need to eventually come up with a more elegant and performant way of handling this type of linking functionality using cloud resources somewhere as it is a fairly common task. With intersections discovered, we can bring everything together into one combined dataset that contains the essential information we will use to create the WikiData items. The following codeblock retrieves the intersection data from PostGIS and sets up a simple function to retrieve the WikiData ID values for those admin entities that I previously assembled whose geometry intersects with ecoregion boundaries.<jupyter_code>df_intersections = get_intersections()
df_ecoregions = recombine_ecoregions()<jupyter_output><empty_output><jupyter_text>In this step, I bring the properties (minus geometry) together again for all ecoregion units across all 5 datasets and apply a couple of processes to prep the data for working up into WikiData items. These include splitting the list of other ecoregion identifiers that will become "part of" relationships in WikiData into a list object and applying our function to provide a list of intersecting administrative boundaries. At this point, we have everything we need to stub out functional WikiData items for each ecoregion unit in all 3 North American classification systems and the 2 for CONUS. This includes linkages we will be able to create to state/province administrative units across all three countries (and each of the countries) along with US counties. We have a representational point location for basic map orientation. And we will be able to build out "part of/has part" relationships between relevant ecoregion units once we have the items created and IDs generated in WikiData. I save out the final processed dataset as a Python data structure for processing to WikiData in another step. At this point, we no longer need our PostGIS intermediary.<jupyter_code>df_ecoregions.head()
cache_final_output(df_ecoregions)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07303238005199716, 'fraction_numerical': 0.015835499881824627, 'mean_word_length': 4.570300157977883, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 5, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '593591', 'license': 'Unlicense', 'n_tokens_mistral': 6665, 'n_tokens_neox': 6361, 'n_words': 2231, 'path': 'ecoregions/Ecoregions Data Processing.ipynb', 'repo_name': 'skybristol/wdbots', 'size': '21155'}