<jupyter_start><jupyter_text>Geo-Spatial processing of GIS shapefiles with Amazon AthenaProcessing geo-spatial data often requires to assign spatial data points to regions or zones, such as zip codes, urban neighborhoods, operational districts, etc.The Athena Engine 2.0 no supports spatial queries.We demonstrate using Athena to join regions or zones to geo-spatial datapoints. The regons tabels are created from ArcGIS shape files.  Data SourcesDownload data from these public data sources: the first set includes geo-spatial locations with longitude and latitude, the remaining three sets are shape files that define areas on the map. We use those to map locations to areas, such as neighborhoods, designated service areas, and zip codes.|Raw data | Description | Source ||---------|-------------|--------|| dockless-vehicles-3.csv | A public dataset of dock-less vehicle rentals, provided by the Office of Civic Innovation and Technology (https://louisvilleky.gov/government/civic-innovation-and-technology/civic-innovation) of the Louisville (KY) Metro Government. The table includes trip data with longitude and latitude of the start and end locations | https://data.louisvilleky.gov/dataset/dockless-vehicles || Dockless_Vehicle_Service_Area.* | Shape files of service areas |  https://data.louisvilleky.gov/dataset/dockless-vehicles || cb_2018_us_zcta510_500k.* | [ZIP Code Tabulation Areas](https://www.census.gov/programs-surveys/geography/guidance/geo-areas/zctas.html) (ZCTAs) are generalized areal representations of United States Postal Service (USPS) ZIP Code service areas. The USPS ZIP Codes identify the individual post office or metropolitan area delivery station associated with mailing addresses. USPS ZIP Codes are not areal features but a collection of mail delivery routes. | https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_zcta510_500k.zip || Louisville_KY_Urban_Neighborhoods | The Louisville Neighborhoods layer consists of polygons representing approximate boundaries and extents of historical and cultural neighborhoods within the pre-merger limits of the City of Louisville, or post-merger Urban Service District. | https://data.lojic.org/datasets/louisville-ky-urban-neighborhoods|   Setup: Load packages, define functionsIn addition to [Pandas](https://pandas.pydata.org) and [NumPy](https://numpy.org) you need to install the following package:1. [s3fs](https://github.com/dask/s3fs/) is required by Pandas. It supports Amazon S3 URIs as file names in the Pandas methords `read_*` and `to_*`. The package loads indirectly with Pandas.2. [PyAthena](https://github.com/laughingman7743/PyAthena/) provides the database connection to Amazon Athena3. [Python Shapefile Library](https://github.com/GeospatialPython/pyshp) reads and writes ESRI Shapefiles.<jupyter_code>import os
import sys
import datetime
import pandas as pd
import numpy as np
import s3fs       # required to read and write tables on Amazon S3. Not called, but importing it ensures that the package is available
import pyathena   # connection to Amazon Athena
import shapefile  # read and write ESRI Shapefiles

S3_BUCKET          = '<Here goes your S3 bucket>'
S3_TABLE_PREFIX    = '<Here goes the S3 prefix>'
S3_ATHENA_QUERY    = f's3://{S3_BUCKET}/{S3_TABLE_PREFIX}/query_output'
ATHENA_SCHEMA_NAME = '<Here goes your Glue database>'
ATHENA_WORK_GROUP  = '<Here goes your Athena Workgroup>'
REGION             = '<Here goes your region>'

def get_athena_connecttion():
    '''
    Open database connection to Amazon Athena
    
    Returns:
        connection object for Pandas
    '''
    return pyathena.connect(s3_staging_dir=S3_ATHENA_QUERY, schema_name=ATHENA_SCHEMA_NAME, region_name=REGION, work_group=ATHENA_WORK_GROUP)
    
    
def save_as_external_athena_table(conn, df, bucket, dirkey, dbname, tblname, overwrite=False):
    '''
    Save Pandas DataFrame as external Athena table. The function supports the following column types: timestamp, int, double, and string.
    The column names will be converted to lower case. The names must not include any special characters, exctept '_'.
    
    Parametes:
        conn:    PyAthena connection
        df:      pandas.DataFrame
        bucket:  S3 bucket
        dirkey:  S3 prefix to data file on S3
        dbname:  Name of GLUE database
        tblname: Name of GLUE table
    '''
    assert dirkey[-1] != '/', "Don't use '/' in S3 key"
    
    df.rename(lambda s: s.lower(), axis=1) \
      .to_csv(f"s3://{bucket}/{dirkey}/table.tsv", index=None, quoting=None, sep='\t', header=None, date_format='%Y-%m-%d %H:%M:%S')
    
    type_renaming = {
        'datetime64[ns]' : 'timestamp',
        'int64' : 'int',
        'float64': 'double',
        'object': 'string',
    }

    cols = [ "`{}` {}".format(r[1], type_renaming[str(r[2])]) for r in df.dtypes.reset_index().itertuples() ]
    cols_string = ',\n'.join(cols)
    
    curs = conn.cursor()
    if overwrite:
        try:
            res = curs.execute(f"""DROP TABLE {dbname}.{tblname}""")
        except pyathena.error.OperationalError as ex:
            if 'Table not found' in str(ex):
                pass
            else:
                raise ex
        
    q = f"""
        CREATE EXTERNAL TABLE IF NOT EXISTS {dbname}.{tblname} (
          {cols_string}
        )
        ROW FORMAT DELIMITED 
        FIELDS TERMINATED BY '\\t' 
        LINES TERMINATED BY '\\n' 
        LOCATION 's3://{bucket}/{dirkey}'
        TBLPROPERTIES ('has_encrypted_data'='false')
    """
  
    res = curs.execute(q)
    return res


def centroid(points):
    '''
    Computes centroid of a 2-dimensional area https://en.wikipedia.org/wiki/Centroid
    
    Parameters:
        points: numpy array of polynom coordinates
        
    Returns:
        x/y-coordinates of centroid, area of centroid
    '''
    pts = np.array(points)
    assert pts.shape[1] == 2, "Not 2-dimensional"
    n  = pts.shape[0]
    x = pts[:,0]
    y = pts[:,1]
    xyxy = [  x[i] * y[(i+1)%n] - x[(i+1)%n] * y[i] for i in range(n) ]
    A = 0.5 * np.sum(xyxy)
    Cx = np.sum([  (x[i] + x[(i+1)%n]) * xyxy[i]  for i in range(n) ])/(6*A)
    Cy = np.sum([  (y[i] + y[(i+1)%n]) * xyxy[i]  for i in range(n) ])/(6*A)
    return Cx, Cy, np.abs(A)


def load_shape_file(shape_file_name):
    '''
    Function to load 
    
    Parameters:
    
    Returns:
        Pandas DataFrame: the field `shape` holds the polygon defintion of the shape that can be interpreted
        by the Athena function `ST_GeometryFromText()`
        Other fields include coordinates of the bounding boxes and centroids, and meta data from the shapefile
    '''
    sf = shapefile.Reader(shape_file_name)
    recs = sf.records()
    shps = sf.shapes()
    fields = sf.fields
    
    shape_df = pd.DataFrame()
    for i, rec in enumerate(recs):
        rdata = {}
        for fld, dt, le, ze in sf.fields[1:]:
            rdata[fld] = [ rec[fld] ]
        bbox = shps[i].bbox
        for j, b in enumerate(['bb_west', 'bb_south', 'bb_east', 'bb_north']):
            rdata[b] = [ bbox[j] ]
        rdata['shape'] = [ """POLYGON(({poly}))""" \
                                .format(poly = ','.join(["{lo} {la}".format(lo=p[0], la=p[1]) 
                                                                        for p in shps[i].points ]))
                         ]
        cog_long, cog_lat, _ = centroid(shps[i].points)
        rdata['cog_longitude'] = cog_long
        rdata['cog_latitude'] = cog_lat

        tmp_df = pd.DataFrame(rdata)
        shape_df = pd.concat([shape_df, tmp_df])
    
    shape_df.index = range(shape_df.shape[0])
    return shape_df<jupyter_output><empty_output><jupyter_text>Process Shape FilesDownload the shape files. Each data set consists of four files with the same name and different extension. Just provide the name of the `.shp` file to load. The other files must be in the same directory. - load shapefiles into Pandas DataFrames with `load_shape_file()`, then- save on S3 and create external table in Athena with `save_as_external_athena_table()`  Zip Codes<jupyter_code>%%time
zip_code_df = load_shape_file('./data/cb_2018_us_zcta510_500k.shp')
print(f"Number of records: {zip_code_df.shape[0]:,}")
display(zip_code_df.head())
%%time
if not 'conn' in globals():
    print("Connecting to Athena")
    conn = get_athena_connecttion()

save_as_external_athena_table(conn, zip_code_df,
                              S3_BUCKET, f"{S3_TABLE_PREFIX}/zipcode_shapes",
                              ATHENA_SCHEMA_NAME, 'zipcode_shapes',
                              overwrite=True)
pd.read_sql("""SELECT * FROM blog1.zipcode_shapes LIMIT 5""", conn)<jupyter_output><empty_output><jupyter_text>Loisville/KY Neighborhoods<jupyter_code>%%time
neighnorhood_df = load_shape_file('./data/Louisville_KY_Urban_Neighborhoods.shp')
print(f"Number of records: {neighnorhood_df.shape[0]:,}")
display(neighnorhood_df.head())
if not 'conn' in globals():
    print("Connecting to Athena")
    conn = get_athena_connecttion()

save_as_external_athena_table(conn, neighnorhood_df,
                              S3_BUCKET, f"{S3_TABLE_PREFIX}/loisville_ky_neighborhoods",
                              ATHENA_SCHEMA_NAME, 'loisville_ky_neighborhoods',
                              overwrite=True)
pd.read_sql("""SELECT * FROM blog1.loisville_ky_neighborhoods LIMIT 5""", conn)<jupyter_output><empty_output><jupyter_text>Dockless Service Zones<jupyter_code>%%time
dockless_zone_df = load_shape_file('./data/Dockless_Vehicle_Distribution_Zones.shp')
print(f"Number of records: {dockless_zone_df.shape[0]:,}")
display(dockless_zone_df)
if not 'conn' in globals():
    print("Connecting to Athena")
    conn = get_athena_connecttion()

save_as_external_athena_table(conn, dockless_zone_df,
                              S3_BUCKET, f"{S3_TABLE_PREFIX}/distribution_zones",
                              ATHENA_SCHEMA_NAME, 'distribution_zones',
                              overwrite=False)
pd.read_sql("""SELECT * FROM blog1.distribution_zones LIMIT 5""", conn)<jupyter_output><empty_output><jupyter_text>Load trip dataUpload the trip data<jupyter_code>trips  = pd.read_csv('data/dockless-vehicles-3.csv')
print(f"Number of records: {trips.shape[0]:,}")
trips.head()
if not 'conn' in globals():
    print("Connecting to Athena")
    conn = get_athena_connecttion()

save_as_external_athena_table(conn, dockless_zone_df,
                              S3_BUCKET, f"{S3_TABLE_PREFIX}/dockless_vehicles",
                              ATHENA_SCHEMA_NAME, 'dockless_vehicles',
                              overwrite=False)<jupyter_output><empty_output><jupyter_text>Athena QueryNow you can try join the neighborhoods to the geo-locations. Learn more about geo-spatial functions in Athena at https://docs.aws.amazon.com/athena/latest/ug/geospatial-functions-list.html<jupyter_code>if not 'conn' in globals():
    print("Connecting to Athena")
    conn = get_athena_connecttion()


q = """
SELECT tr.startdate, tr.starttime, tr.endtime, tr.tripduration, tripdistance 
      , nb1.nh_code AS start_nbid, nb1.nh_name AS start_neighborhood
      , nb2.nh_code AS end_nbid, nb2.nh_name AS end_neighborhood
  
FROM "{database}"."dockless_vehicles" tr

JOIN "{database}"."loisville_ky_neighborhoods" nb1
      ON ST_Within(ST_POINT(CAST(tr.startlongitude AS DOUBLE), CAST(tr.startlatitude AS DOUBLE)), ST_GeometryFromText(nb1.shape))  
JOIN "{database}"."loisville_ky_neighborhoods" nb2
      ON ST_Within(ST_POINT(CAST(tr.endlongitude AS DOUBLE), CAST(tr.endlatitude AS DOUBLE)), ST_GeometryFromText(nb2.shape))

LIMIT 100;
""".format(
    database=ATHENA_SCHEMA_NAME
)

trips_with_nbhd = pd.read_sql(q, con=conn)
print(f"Number of records: {trips_with_nbhd.shape[0]:,}")
trips_with_nbhd.head()<jupyter_output>Number of records: 100
{'fraction_non_alphanumeric': 0.10978398983481576, 'fraction_numerical': 0.011012282930961457, 'mean_word_length': 3.3119065010956903, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 13, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '604977', 'license': 'MIT-0', 'n_tokens_mistral': 4027, 'n_tokens_neox': 3815, 'n_words': 1004, 'path': 'notebook/GeoSpatialProcessingGISshapefilesWithAmazonAthena.ipynb', 'repo_name': 'aws-samples/dynamodb-ml-prediction-amazon-athena', 'size': '11805'}