# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Find pools in City of Redlands

# +
import os
from arcgis.gis import GIS
from arcgis.geocoding import geocode

gis = GIS('https://deldev.maps.arcgis.com', 'demos_deldev', os.environ.get('DEMOSPWD'))
# -

# # Inputs

extent = geocode('City of Redlands, CA', out_sr=3857)[0]['extent']
extent

zoning_lyr = gis.content.get('d83b7a006f444212ac81a16399a37ac6') # Redlands Zoning 
zoning_lyr

naip = gis.content.get('e4da3b6720f545aeaaf3fe8141da1e21') # USA NAIP Imagery: Color Infrared
naip

# # Deploy trained model

# + [markdown] heading_collapsed=true
# ### Load model

# + hidden=true
import torch
import torch.nn as nn
from PIL import Image
import numpy as np

from util import load_model, norm, denorm, export_img, get_tile_images, get_img, \
        get_cropped_tiles, predict_, detect_objects, suppress_close_pools, overlap, predict_classf
from model import ConvnetBuilder, SSD_MultiHead, resnet34, k

torch.backends.cudnn.benchmark=True

f_model = resnet34
pretrained = '9.h5'

# + hidden=true
head_reg = SSD_MultiHead(k, -4.)
models = ConvnetBuilder(f_model, 0, 0, 0, custom_head=head_reg)
model = models.model
model = load_model(model, pretrained)
model.eval()

path = r'inferencing/'

# + hidden=true
# !mkdir $path

# + hidden=true
import os
import math
import pandas as pd

from arcgis.raster.functions import apply
from arcgis.features import SpatialDataFrame
from arcgis.mapping import MapImageLayer
from arcgis.geometry import Envelope, project
from IPython.display import display

ext = Envelope(extent)
city_ext = project([ext], 3857, 3857)[0]

# + hidden=true
naiplayer = naip.layers[0]
naipfalse = apply(naiplayer, 'FalseColorComposite')

# + hidden=true
df = SpatialDataFrame.from_layer(zoning_lyr.layers[0])
residential_areas = list(df[df.Zoning=='Residential']['SHAPE'].extent)

areas = [Envelope({"xmin": x[0], "ymin": x[1], "xmax": x[2], "ymax": x[3], "spatialReference": 3857}) for x in residential_areas]

# + hidden=true
sz = 224
pixx, pixy = 90, 90
N = 16
w = h = N*pixx - pixx/2

nX = math.ceil((city_ext['xmax'] - city_ext['xmin'])/w)
nY = math.ceil((city_ext['ymax'] - city_ext['ymin'])/h)

years      = [2014, 2016]
objects    = []
inferenced = []
done_imgs  = {}

borders = [ (45,   45,   45,   45),   # xmin, ymin, xmax, ymax - map space
          (22.5, 22.5, 67.5, 67.5),   # bottom left corner
          (22.5, 67.5, 67.5, 22.5),   # top left corner
          (67.5, 22.5, 22.5, 67.5),   # bottom right corner
          (67.5, 67.5, 22.5, 22.5) ]  # top right corner

crops   = [ (112, 112, 112, 112),  # left, top, right, bottom -> xmin, ymax, xmax, ymin
            ( 56, 168, 168,  56),  # bottom left corner
            ( 56,  56, 168, 168),  # top left corner
            (168, 168,  56,  56),  # bottom right corner
            (168,  56,  56, 168) ] # top right corner
# -

# ### Detect swimming pools

for year in years:
    for i in range(nX):
        for j in range(nY):
            filename = f'redlands_large_{i}_{j}_({year}).jpg'

            xmin = city_ext['xmin'] + (j)*w
            ymin = city_ext['ymax'] - (i+1)*h
            xmax = xmin + N*pixx
            ymax = ymin + N*pixy

            ext = {'xmin':xmin, 'ymin':ymin, 'xmax':xmax, 'ymax':ymax, 'spatialReference':3857}
            thisarea = Envelope(ext)

            for area in areas:
                if (overlap(area, ext)):
                    print(filename, 'Detecting...')

                    imgpath = path + filename

                    if not os.path.isfile(imgpath):
                        print(filename, 'Downloading...')
                        imgpath = naipfalse.export_image(ext, image_sr={'wkid': 3857}, 
                                         bbox_sr=3857, size=[224*8, 224*8], f='image', 
                                         export_format='jpg', adjust_aspect_ratio=False,
                                         mosaic_rule={"mosaicMethod": "esriMosaicAttribute", 
                                                "where": f'Year = {year}', 
                                                "sortField":"Year", 
                                                "sortValue":"3000", 
                                                "ascending":True, 
                                                "mosaicOperation":"MT_FIRST"},
                              save_folder=path,
                              save_file=filename)
                    
                    if filename not in inferenced:                        
                        img = np.array(Image.open(imgpath))

                        tiles = get_tile_images(img)
                        pools = detect_objects(model, tiles, ext, cycle=0)
                        print(filename, 'Found ' + str(len(pools)) + ' pools (pass #0)')
                        objects.extend(pools)

                        for iii, crop in enumerate(crops):
                            paddedtiles = get_cropped_tiles(img, crop)
                            pools = detect_objects(model, paddedtiles, ext, crop=borders[iii], cycle=iii+1)
                            print(filename, 'Found ' + str(len(pools)) + f' pools (pass #{iii+1})')
                            objects.extend(pools)

                        inferenced.append(filename)

# # Post processing

# + [markdown] heading_collapsed=true
# ### Create feature layers

# + hidden=true
import pickle 
pickle.dump(objects, open(path + 'objects.pkl', 'wb'))

# + hidden=true
import pickle
objects = pickle.load(open(path + 'objects.pkl', 'rb'))

# + hidden=true
len(objects)

# + hidden=true
for i, obj in enumerate(objects):
    obj['SHAPE'] = {'x': obj['x'], 'y': obj['y']}

sdf = SpatialDataFrame.from_dict(objects)
sdf = sdf.reset_index()
sdf.sr = {'wkid': 3857}

# + hidden=true
selected, suppressed = suppress_close_pools(sdf, 15)

# + hidden=true
selected.to_featureclass(r'inferencing', 'selected_pools.shp')

# + hidden=true
selected.to_featureclass(r'inferencing', 'suppressed_pools.shp')

# + hidden=true
import matplotlib.pyplot as plt
import numpy as np

# %matplotlib inline
plt.plot([len(selected[selected.score > i]) for i in np.arange(0, 100)]);

# + hidden=true
print('No. of pools detected: ', len(selected[selected.score > 30]))

# + hidden=true
gis = GIS('https://deldev.maps.arcgis.com', 'demos_deldev', os.environ.get('DEMOPWD'))

# + hidden=true
pools = gis.content.import_data(selected, title='Pools in Redlands', tags='deeplearning')    
pools

# + hidden=true
aoi_parcels = parcels.layers[0]
# pools = gis.content.get('6531c1e0d3234bb8a7c630302060a7c3')

# + hidden=true
from arcgis.mapping.renderer import generate_renderer

parcels = gis.content.get('fd0ab9419c234a4783da3466e804ce85')
parcels_with_unauthorized_pools = gis.content.get('d02e7a19fdeb46c5ace518535c0c9646')
detected_pools = gis.content.get('dae7b7c72e0c41f88852597ec517fca1')

def show_results(rm):
    red_renderer = generate_renderer('Polygon', symbol_type='simple', symbol_style='s',
                                     colors=[[0,0,0,0]], outline_color=[255,0,0,255], line_width=2.25)
    rm.add_layer(parcels, {"renderer":{"type":"uniqueValue","field1":"pool_","defaultSymbol":{"color":None,"outline":{"color":[153,153,153,255],"width":0.75,"type":"esriSLS","style":"esriSLSSolid"},"type":"esriSFS","style":"esriSFSSolid"},"defaultLabel":"Other","uniqueValueInfos":[{"value":"UNH","symbol":{"color":None,"outline":{"color":[0,255,255,255],"width":1.125,"type":"esriSLS","style":"esriSLSSolid"},"type":"esriSFS","style":"esriSFSSolid"},"label":"UNH"},{"value":"HEA","symbol":{"color":None,"outline":{"color":[0,255,255,255],"width":1.125,"type":"esriSLS","style":"esriSLSSolid"},"type":"esriSFS","style":"esriSFSSolid"},"label":"HEA"}]}})
    rm.add_layer(parcels_with_unauthorized_pools, {'renderer': red_renderer})
    #result_map.add_layer(detected_pools, {'renderer': { "type": "uniqueValue", "field1": None,"valueExpression": "if (($feature.pool_ == 'UNH') || ($feature.pool_ == 'HEA')) {\n    return \"Authorized Pool\";\n} else if ($feature.score > 75) {\n    return \"Unauthorized Pool (high confidence)\";\n} else  if ($feature.score > 60) {\n    return \"Unauthorized Pool (low confidence)\";\n}","valueExpressionTitle": "Custom","uniqueValueInfos": [{"value":"Authorized Pool","symbol":{"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS","url":"https://rohitgeo.github.io/images/ArrowYellowCenter.png","width":31.5,"height":31.5},"label":"Authorized Pool"},{"value":"Unauthorized Pool (high confidence)","symbol":{"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS","url":"https://rohitgeo.github.io/images/ArrowRedCenter.png","width":31.5,"height":31.5},"label":"Unauthorized Pool (high confidence)"},{"value":"Unauthorized Pool (low confidence)","symbol":{"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS","url":"https://rohitgeo.github.io/images/ArrowOrangeCenter.png","width":31.5,"height":31.5},"label":"Unauthorized Pool (low confidence)"}]}})
    
arcade_renderer = {'renderer': { "type": "uniqueValue", "field1": None,"valueExpression": "if (($feature.pool_ == 'UNH') || ($feature.pool_ == 'HEA')) {\n    return \"Authorized Pool\";\n} else if ($feature.score > 75) {\n    return \"Unauthorized Pool (high confidence)\";\n} else  if ($feature.score > 60) {\n    return \"Unauthorized Pool (low confidence)\";\n}","valueExpressionTitle": "Custom","uniqueValueInfos": [{"value":"Authorized Pool","symbol":{"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS","url":"https://rohitgeo.github.io/images/ArrowYellowCenter.png","width":31.5,"height":31.5},"label":"Authorized Pool"},{"value":"Unauthorized Pool (high confidence)","symbol":{"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS","url":"https://rohitgeo.github.io/images/ArrowRedCenter.png","width":31.5,"height":31.5},"label":"Unauthorized Pool (high confidence)"},{"value":"Unauthorized Pool (low confidence)","symbol":{"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS","url":"https://rohitgeo.github.io/images/ArrowOrangeCenter.png","width":31.5,"height":31.5},"label":"Unauthorized Pool (low confidence)"}]}}
   

# + [markdown] heading_collapsed=true
# ### Analysis: Join pools to residential parcels

# + hidden=true
parcels = gis.content.get('fd0ab9419c234a4783da3466e804ce85') # Parcels_In_Redlands 
parcels

# + hidden=true
from arcgis.features.summarize_data import join_features

pools_within_parcels = join_features(pools, aoi_parcels,
                                           join_operation='JoinOneToMany',
                                           spatial_relationship='intersects',
                                           output_name="Pools_Within_Parcels")
pools_within_parcels

# + [markdown] hidden=true
# #### Find pool with highest score in parcel

# + hidden=true
df = pools_within_parcels.layers[0].query().df
df.to_featureclass(r'inferencing', 'Pools_Within_Parcels.shp')

# + hidden=true
df = SpatialDataFrame.from_featureclass('inferencing/Pools_Within_Parcels.shp')

# + hidden=true
grouped = (df[(df['score'] > 30) & (df['c_score'] > 50.0)].groupby(by='apn').apply(lambda x: x.loc[x.c_score.idxmax(), :]))
grouped.reset_index(drop=True, inplace=True)

temp = SpatialDataFrame(data=grouped)
SHAPE = temp['SHAPE'].copy()
del temp['SHAPE']
temp = temp.set_geometry(SHAPE)
temp.sr = {'wkid': 3857}

detected_pools = gis.content.import_data(temp, title='Detected Swimming Pools', tags='deeplearning')    
detected_pools

# + hidden=true
unauthorized_pools = detected_pools.layers[0]
unauthorized_pools.filter = "(pool_ <> 'HEA') AND (pool_ <> 'UNH') AND (score > 60)"
# -

# ### Analysis: Find parcels with unassessed swimming pools

parcels_with_unassessed_pools = join_features(aoi_parcels, unassessed_pools, 
                                    join_operation='JoinOneToMany',
                                    spatial_relationship='intersects',
                                    output_name="Parcels_With_Unassessed_Pools")
parcels_with_unassessed_pools

# + [markdown] heading_collapsed=true
# ### Create web map

# + hidden=true
result_map = gis.map('Redlands County Club, Redlands, CA', zoomlevel=17)
result_map.basemap = 'satellite'
result_map

# + hidden=true
show_results(result_map)

# + hidden=true
webmap_item = result_map.save({
    'title': 'Detected swimming pools',
    'tags': 'deeplearning',
    'snippet': 'Swimming pools detected using deep learning'
})
webmap_item

# + [markdown] heading_collapsed=true
# # Generate Report

# + hidden=true
parcel_df = parcels_with_unassessed_pools.layers[0].query().df

# + hidden=true
parcel_df[['apn', 'pool_', 'spa_', 'taxstatus', 'landvalue', 'improvemen']] 

# + [markdown] hidden=true
# ## Looking in closer

# + hidden=true
m = gis.map('1418 Sterling Rd, Redlands, CA')
m.basemap = 'satellite'
m

# + hidden=true
show_results(m)
m.add_layer(detected_pools, arcade_renderer)
# -

# # Zombie pools

# + [markdown] heading_collapsed=true
# ## Classify pools as 'clean' or 'green' using deep learning

# + hidden=true
folder = 'inferencing/cleangreen_2706'

# + hidden=true
imagery = gis.content.get('10df2279f9684e4a9f6a7f08febac2a9')
imagery

# + hidden=true
imagery_lyr = MapImageLayer.fromitem(imagery)

# + hidden=true
for index, row in df.iterrows():
    img_name = 'img' + str(df.loc[index,'parcelnumb']).zfill(4)# index).zfill(4)
    try:
        g = row['SHAPE']
        x, y = g['x'] - 7, g['y'] - 7
        ext = (x, y, x + 14, y + 14)
        if not os.path.isfile(folder + '/' + img_name + '.jpg'):
            imagery_lyr.export_map(ext, 102100, '224,224', f='image', format='jpg',
                    save_folder=folder, save_file='{}.jpg'.format(img_name))
            print(f'{img_name} downloaded..')    
    except Exception as e:
        print('Error {} getting {}'.format(e,img_name))

# + [markdown] heading_collapsed=true hidden=true
# ### Load model for classifying zombie pools

# + hidden=true
model_path = 'clean_green.h5'

# + hidden=true
models = ConvnetBuilder(resnet34, 2, False, False, pretrained=True)
clean_green_model = models.model
clean_green_model = load_model(clean_green_model, model_path)

# + hidden=true
clean_green_model.eval();

# + hidden=true
cleangreen_scores = {}

for index, row in df.iterrows():
    
    try:
        img_name = 'img' + str(index).zfill(4) + '.jpg'
        img_path = folder + '/' + img_name
        img = np.array(Image.open(img_path))
        score = float(predict_classf(clean_green_model,
                               norm(img).transpose(2, 0, 1))[0])
        print(img_name, score)
        cleangreen_scores[index] = score
        
    except Exception as e:
        cleangreen_scores[index] = -1
        print(img_name, '-1')

# + hidden=true
df['neglected'] = pd.Series(list(cleangreen_scores.values())) * 100.0

# + hidden=true
df

# + hidden=true
df.sr = {'wkid': 3857}

# + [markdown] heading_collapsed=true hidden=true
# ### Create layer of pools needing inspection

# + hidden=true
zombiepools_lyr = portal.content.import_data(df, title='Pools_needing_maintenance_in_Redlands (27 Jun 2018-a)', tags='datascience')    
zombiepools_lyr

# + hidden=true
inspection_pools = df[(df.neglected > 76.5) & (df.score > 45)]

# + hidden=true
inspection_lyr = portal.content.import_data(inspection_pools, title='Pools_needing_inspection (27 Jun 2018-a)', tags='datascience')    
inspection_lyr
# -

# ## Create layer of pools needing inspection

# +
from arcgis.features.summarize_data import join_features

parcels_needing_inspection = join_features(aoi_parcels, inspection_lyr,
                                    join_operation='JoinOneToMany',
                                    spatial_relationship='intersects',
                                    output_name="Parcels_needing_inspection")
parcels_needing_inspection

# + [markdown] heading_collapsed=true
# # Imagery Visit app to inspect results

# + hidden=true
from IPython.display import IFrame
IFrame('http://maps.arcgis.com/apps/ImageVisit/index.html?appid=5502f61d9d944cd2900021970682bfed', width=900, height=500)

# + [markdown] heading_collapsed=true
# # Create assignments for Workforce for ArcGIS

# + hidden=true
pools_layer = gis.content.get("571f435be7dc43b3aaf24aabec96070a").layers[0]

# get the pools that require inspections
pools = pools_layer.query("score > 70.7 AND c_score > 200")

# + hidden=true
# The Workforce project
project = Project(gis.content.get("9ee1b2a4cb884ecaaa7ef350a77db04d"))
# The first/only dispatcher (demo_deldev)
dispatcher = project.dispatchers.get(1)
# The assignment type
pool_inspection = project.assignment_types.get(name="Pool Inspection")

# list of assignments to add
assignments = []
for pool in pools:
    
    # Set priority based on neglected field (default to "none")
    priority = "none"
    if pool.attributes["neglected"] >= 90:
        priority = "critical"
    elif pool.attributes["neglected"] >= 80:
        priority = "high"
    elif pool.attributes["neglected"] >= 70:
        priority = "medium"
    elif pool.attributes["neglected"] >= 60:
        priority = "low"
    
    # Create the assignment objects to add
    assignments.append(
        Assignment(
            project,
            geometry=pool.geometry,
            location=pool.attributes["ownername"], # This is a string that represent the location/address
            status="unassigned",
            description="Pool Inspection Required", # Any description (could build a string based on other fields)
            assignment_type=pool_inspection,
            dispatcher=dispatcher,
            priority=priority,
            due_date=datetime.now()+timedelta(days=2) # This attribute is not required and can be omitted if desired
        )
    )
    
# Add the all of the assignments
assignments = project.assignments.batch_add(assignments)

StarcoderdataJupyterScriptsDedupFiltered