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

# + [markdown] id="lYeNN6UJDMDz" colab_type="text"
# # **Population prediction via CNNs and geospatial imagery** #
# ### Spring 2020 | AM216 Final Project
# ### <NAME>, <NAME>, <NAME>
# GitHub repository [here](https://github.com/chickert/geospatial_analysis)
#

# + [markdown] id="mgtEVoQAc5BL" colab_type="text"
# #Abstract:
#
# Google Earth Engine has a multi-petabyte catalog of open source satellite imagery and geospatial datasets available to scientists and researchers to analyze. One interesting problem is extracting demographic information from these satellite images, such as the population of a given area. Our goal is to build neural network models to estimate the population from satellite images. This is a difficult task, as the model must learn land use types within an image in order to identify high- versus low-population areas. Previous sources have worked on a similar problem using deep sequential models and UNet architectures. We expand this by building many models including FFNNs, CNNs, UNets, and also investigating the value of transfer learning for this task.

# + [markdown] id="lNyso3y1EMAF" colab_type="text"
# # Introduction #
# Our group is interested in the combination of machine learning methods and satellite imagery. The amount of data in this space (pun intended) is already remarkable, and growing rapidly. In addition to the large datasets released from government-affiliated science organizations -- such as NASA/USGS' Landsat image collection or the EU Copernicus program's Sentinel data hub -- the proliferation of privately-held "constellations" is now contributing to the growing data pool. Simultaneously, improvements in image analysis driven by machine learning methods make gleaning insights from this relatively-unstructured data more feasible than ever. 
#
# For our project, our group assesses the performance of neural networks in predicting a US census tract's population with only raw RGB pixel data collected from the Copernicus satellite program's geospatial imagery. For the context of this project, we will limit our anlaysis to Massachusetts, though we expect that our results are relevant beyond the state's borders. On an even broader level, we hope that the general process and data pipeline used to implement our project could be adapted relatively easily to a much wider variety of predicton tasks using geospatial data from Google's Earth Engine. 
#

# + [markdown] id="5ci63qNtFHi8" colab_type="text"
# # Background and Data Sources #
# ## 2010 Census Data ##
# Every 10 years, the United States conducts a census and releases this information to the public. The highest level of granularity in this information is a "census block", which is roughly the size of a city block. Census blocks are, in turn, part of census tracts, which vary significantly in size depending on population density, but in urban areas, are roughly the size of a neighborhood and contain somewhere between [1,200 and 8,000 inhabitants](https://www2.census.gov/geo/pdfs/education/CensusTracts.pdf).

# + [markdown] id="h1wKL6lTIUYf" colab_type="text"
# ## Cenpy ##
#
# The 2010 Census is massive, and as such, we used CenPy (https://cenpy-devs.github.io/) to easily query this dataset for the information we wanted. CenPy is built with Geopandas (https://geopandas.org/), which allows for easy data manipulation and discovery.
#
# Below, we install CenPy and pull back information about every tract in the state of Massachusetts. 
#
# Specifically, note the variable paramater, which controls which data we get back from the census, for example `variables=['H006002', 'P001001']`. Every variable is defined [here](https://api.census.gov/data/2010/dec/sf1/variables.html). In this example the code is pulling back
# `H006002` which corresponds to "Average household size!!Renter occupied"
# and `P001001` which corresponds to "Total Population". In the code block below we pull back the population data using this second variable.
#
# It should again be noted that aggregating this information using the raw census data would be slow and somewhat computationally expensive; thus, using CenPy is preferred.
#

# + id="zZRKoxcdDEf2" colab_type="code" outputId="78f0520c-9f75-4a52-f17e-118c0f7ddce7" colab={"base_uri": "https://localhost:8080/", "height": 928}
# !apt install libspatialindex-dev
# !pip install cenpy
import cenpy
massachusetts = cenpy.products.Decennial2010().from_state('Massachusetts',  variables=['P001001'], level="tract")
massachusetts.head()


# + [markdown] id="-CiHaoXTELUJ" colab_type="text"
# Now we filter this to only look at one of our group member's Census tract in Cambridge, MA, which is in `state = 25`, `county=017`, `tract = 353300`
#
# ---
#
#
#

# + id="YrkDR_gpMUlD" colab_type="code" outputId="103cea0e-9249-4b83-91de-e8c9269b64e2" colab={"base_uri": "https://localhost:8080/", "height": 94}
massachusetts[(massachusetts['state'] == '25') & (massachusetts['county'] == '017') & (massachusetts['tract'] == '353300')]

# + [markdown] id="jFbKFhtkDDmn" colab_type="text"
# Thus we see that this census tract has a total population of 3,471. 

# + [markdown] id="6vzuTiWUNVj0" colab_type="text"
# ## Google Earth Engine
# Google Earth Engine allows for easy access to a vast array of satellite imagery, and -- importantly for our project -- also contains geometries of Census Tracts and Census Blocks. Below we authenticate Earth Engine to run within this Colab notebook.
#

# + id="Lz9Sa5_cNRaY" colab_type="code" colab={}
import ee

# + id="ZjdpPkryONSk" colab_type="code" outputId="e47413d9-9282-487a-be10-dc1e429f6411" colab={"base_uri": "https://localhost:8080/", "height": 178}
# Trigger the authentication flow.
ee.Authenticate()

# Initialize the library.
ee.Initialize()

# + [markdown] id="2-RJcFstPCko" colab_type="text"
# We will also allow this to read and write from Google Drive, and then import the `folium` library to visualize maps.

# + id="yMZ_XekmOP_W" colab_type="code" colab={}
from google.colab import auth
auth.authenticate_user()

# + id="4uFQQV4iOqIu" colab_type="code" outputId="9257938f-7e30-4834-f8c4-3e0a52c881b6" colab={"base_uri": "https://localhost:8080/", "height": 34}
import folium
print(folium.__version__)

# + [markdown] id="b0s-VNmaV59n" colab_type="text"
# ## Satellite Images
#

# + [markdown] id="-pniJfM1QlIZ" colab_type="text"
# ### Landsat
# The USGS and NASA's Landsat satelites have been taking pictures of the world since the 1970s. The two most recent satellites, Landsat 7 and 8, are both integrated with Earth Engine, and provide a vast array of imformation such a surface reflection, burn, vegatation, and temperature information.
#
# Ultimately -- and as we describe below -- we elected to go with data from the EU's Copernicus program. However, it is instructive to first provide some insight into our investigation of the Landsat data (especially for readers considering alternative projects with geospatial data), so we provide that description below. See further below for the data processing we do with the Copernicus dataset, which we ultimately used for our project. 
#
# In both cases, we rely strictly on the RGB bands from the surface reflection data, which correspond to bands 2,3, and 4.
#
# Given our Census data is from 2010, we initially wanted to get data from as close to 2010 as possible, which would come from Landsat 7, which has been recording images since 1999.
#

# + [markdown] id="j44t9B5iV-Tq" colab_type="text"
#
# ### Removing Clouds
# One consideration when working with satellite images is the presense of clouds in many images. To combat this, we mask them via bitshifting with a method  adaped from this [notebook](https://colab.research.google.com/github/google/earthengine-api/blob/master/python/examples/ipynb/AI_platform_demo.ipynb).
#

# + id="tu8QG_lLP7Hi" colab_type="code" outputId="e3ff23df-21a8-41df-b284-151281bdbce1" colab={"base_uri": "https://localhost:8080/", "height": 674}
# Cloud masking function.
fc = ee.FeatureCollection('TIGER/2016/States').filter(ee.Filter.eq('NAME', 'Massachusetts'))
test = fc.getInfo()['features']
feature = ee.Feature(test[0])
aoi = feature.geometry()
L8SR = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
BANDS = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']

# Cloud masking function.
# Cloud masking function.
def maskL8sr(image):
  cloudShadowBitMask = ee.Number(2).pow(3).int()
  cloudsBitMask = ee.Number(2).pow(5).int()
  qa = image.select('pixel_qa')
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(
    qa.bitwiseAnd(cloudsBitMask).eq(0))
  return image.updateMask(mask).select(BANDS).divide(10000)

# The image input data is a 2018 cloud-masked median composite.
image = L8SR.filterDate('2013-01-01', '2013-12-31').map(maskL8sr).median().clip(aoi)

# Use folium to visualize the imagery.
mapid = image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3})
map = folium.Map(location=[42.3612549,-71.1088637])

folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='median composite',
  ).add_to(map)
map.add_child(folium.LayerControl())
map

# + [markdown] id="43emI_KzVUDZ" colab_type="text"
# Note the map above, which has a map of Massachusetts with clouds removed.
#
# Now we will examine the specific census tract in Cambridge, MA that we pinpointed before. The geometries of every US census tract are accessible via https://developers.google.com/earth-engine/datasets/catalog/TIGER_2010_Tracts_DP1, which also includes some additional information about each tract.

# + id="lPYIp3GmW-G-" colab_type="code" outputId="ad6c632a-2109-462c-a7b9-e6c7625eaf9d" colab={"base_uri": "https://localhost:8080/", "height": 674}

# Cloud masking function.
fc = ee.FeatureCollection('TIGER/2016/States').filter(ee.Filter.eq('NAME', 'Massachusetts'))
test = fc.getInfo()['features']
feature = ee.Feature(test[0])
aoi = feature.geometry()
L8SR = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
BANDS = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']

# Cloud masking function.
# Cloud masking function.
def maskL8sr(image):
  cloudShadowBitMask = ee.Number(2).pow(3).int()
  cloudsBitMask = ee.Number(2).pow(5).int()
  qa = image.select('pixel_qa')
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(
    qa.bitwiseAnd(cloudsBitMask).eq(0))
  return image.updateMask(mask).select(BANDS).divide(10000)

# The image input data is a 2018 cloud-masked median composite.
image = L8SR.filterDate('2013-01-01', '2013-12-31').map(maskL8sr).median().clip(aoi)

# Use folium to visualize the imagery.
mapid = image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3})
map = folium.Map(location=[42.3612549,-71.1088637])

folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='median composite',
  ).add_to(map)
map.add_child(folium.LayerControl())
map

# + id="jEwsKKnsVe_c" colab_type="code" outputId="2e536dbb-4faa-48f3-cc15-a52e9931e1c6" colab={"base_uri": "https://localhost:8080/", "height": 710}
fc = ee.FeatureCollection("TIGER/2010/Tracts_DP1").filter(ee.Filter.eq('geoid10', '25017353300'))
print(fc.size().getInfo())
test = fc.getInfo()['features']
feature = ee.Feature(test[0])
aoi = feature.geometry()
print(fc.size().getInfo())
L8SR = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
BANDS = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']

# Cloud masking function.
# Cloud masking function.
def maskL8sr(image):
  cloudShadowBitMask = ee.Number(2).pow(3).int()
  cloudsBitMask = ee.Number(2).pow(5).int()
  qa = image.select('pixel_qa')
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(
    qa.bitwiseAnd(cloudsBitMask).eq(0))
  return image.updateMask(mask).select(BANDS).divide(10000)

# The image input data is a 2018 cloud-masked median composite.
image = L8SR.filterDate('2013-01-01', '2013-12-31').map(maskL8sr).median().clip(aoi)

# Use folium to visualize the imagery.
mapid = image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.3})
map = folium.Map(location=[42.3612549,-71.1088637], zoom_start=15)

folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='median composite',
  ).add_to(map)
map.add_child(folium.LayerControl())
map


# + [markdown] id="XKhSefdbYQEG" colab_type="text"
# Above is the census tract. Unfortunately, this image is much blurrier than anticipated. A quick investigation reveals that Landsat has a resolution of 30 meters per pixel, which seems suboptimal for the purpose of pulling information from at the lowest level. As such, we decided to move to another satellite image source found on Earth Engine, the EU's Copernicus Sentinel-2 satellite (https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/resolutions/radiometric), which has a resolution of 10 meters per pixel, which we believe is better for our analysis.
#
# It should be noted that the images that Google uses for Google Earth and Google Maps are **not** licensed for use on Earth Engine, and are provided by DigitalGlobe. These images are as good as 0.31 meters per pixel, but are not free for use and are licensed by Google. As we detail below, accessing this data would be one path worth pursuing in future projects to improve the models' performance. 

# + [markdown] id="bQs6SE2YaA6E" colab_type="text"
# ### Copernicus
# Cambridge census tract using the [Copernicus Sentinel-2 Surface Reflection data](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR) as the image source.

# + id="TcoG9IUeXnaG" colab_type="code" outputId="349028b3-6bc0-49ec-b3be-6fdf7f38163b" colab={"base_uri": "https://localhost:8080/", "height": 710}
fc = ee.FeatureCollection("TIGER/2010/Tracts_DP1").filter(ee.Filter.eq('geoid10', '25017353300'))
print(fc.size().getInfo())
test = fc.getInfo()['features']
feature = ee.Feature(test[0])
aoi = feature.geometry()
print(fc.size().getInfo())
image = ee.ImageCollection('COPERNICUS/S2_SR').filter(ee.Filter.calendarRange(4,10,'month')).filterBounds(aoi).filterMetadata('CLOUD_COVERAGE_ASSESSMENT', 'less_than', 30)
image = image.median().clip(aoi)

mapid = image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 2000})
map = folium.Map(location=[42.3612549,-71.1088637], zoom_start=15)

folium.TileLayer(
    tiles=mapid['tile_fetcher'].url_format,
    attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='median composite',
  ).add_to(map)
map.add_child(folium.LayerControl())
map


# + [markdown] id="n3DXFdgcdjk7" colab_type="text"
# As one can see, this image is significantly higher resolution than the Landsat 7 image. Moreover, the process of removing cloud is easier as Copernicus computes a `CLOUD_COVERAGE_ASSESSMENT`. Thus, in the code above, we filter to focus on the census block, and then take the median of all images of that geography with a `CLOUD_COVERAGE_ASSESSMENT` of less than 30. 
#
# One limitation here is that images are not availble from the earlier part of of the 2010s. Thus, one key issue with our subsequent analyses is that that the images are not from the same year as the data.

# + [markdown] id="9Z0ssUMo5mVW" colab_type="text"
# # Data Pipeline
# ## Saving Images
# While Earth Engine does provide native functionality for exporting images to Tensorflow/Keras, we found it was actually easier to simply export each image as a PNG, thus creating a layer abstraction between our models and the data itself. This also made it easier to share the dataset amongst our group members, and the .png files made for a familiar format for use in our models. We exported the images in the following manner, labeling each tract by its `geoid`.
#
#

# + [markdown] id="AQUl5YNY6vEs" colab_type="text"
# First, pull back each tract's data

# + id="L-pdqqVV6st2" colab_type="code" outputId="8ca57862-8c54-41d5-fe9c-499a03fd76ce" colab={"base_uri": "https://localhost:8080/", "height": 893}
# !apt install libspatialindex-dev
# !pip install cenpy
import cenpy
massachusetts = cenpy.products.Decennial2010().from_state('Massachusetts',  variables=['P001001'], level="tract")
massachusetts.head()

# + [markdown] id="w2FbKk-t6_2M" colab_type="text"
# Below we download a 512 x 512 version of the images. 

# + id="vR-VlD88bF4N" colab_type="code" outputId="839d3559-6c9e-4692-8779-6374fa686c37" colab={"base_uri": "https://localhost:8080/", "height": 493, "referenced_widgets": ["53a91230688a4d3988aeb9387e9a9150", "e89b319e273f46e78fa70db34de85fd5", "d486a91fcc9848cc81ecaa17d7bf3054", "c82da68db2634ccfab978502105a1a05", "b417ba7d54c14cb6b4c1da056563dcc6", "64382460351c48419fa3b7104539fe48", "e198b2fa372a4b2c8bce49b85857188f", "b6b96281e5664fe0902c0bc5418ca1b5"]}
import os.path
from tqdm.auto import tqdm
grouped = massachusetts.groupby(['tract', 'county'])
tracts = grouped[['tract', 'county']]
def save_tract_image(tract):
  (this_tract, this_county) = tract[0]
  geoid = '25'+this_county+this_tract
  save_path = "/content/drive/My Drive/summer_retry/"+geoid+".png"
  if os.path.exists(save_path):
    pass
  else:
    fc = ee.FeatureCollection("TIGER/2010/Tracts_DP1").filter(ee.Filter.eq('geoid10', geoid))
    test = fc.getInfo()['features']
    feature = ee.Feature(test[0])
    aoi = feature.geometry()
    image = ee.ImageCollection('COPERNICUS/S2_SR').filter(ee.Filter.calendarRange(4,10,'month')).filterBounds(aoi).filterMetadata('CLOUD_COVERAGE_ASSESSMENT', 'less_than', 30)
    image = image.median().clip(aoi)
    parameters = {'min': 0, 
                'max': 2000,
                'dimensions': 512,
                'bands': ['B4', 'B3', 'B2'],
                'region':aoi}
    url = image.getThumbUrl(parameters)
    # urllib.request.urlretrieve(url, "/content/drive/My Drive/summer_retry/"+geoid+".png")


for tract in tqdm(tracts):
  save_tract_image(tract)

# + [markdown] id="YbFsy05K_vFX" colab_type="text"
# One potential issue with this approach is that we create a 512x512 image of every census tract, despite the fact that tracts' land areas can differ significantly. One thing we can do to potentially combat this is to also record the area of each tract, and later feed that as another parameter into the models. We try this later in the experiments section. Thus, below we will create a dictionary of each tract and and its area, as calculated by the the area of the geometry in Earth Engine (which has a convenient built-in method to do so). We divide by `1000 * 1000` to convert the area from square meters to square kilometers. 

# + id="RQkL1OHi7FDA" colab_type="code" outputId="6d799fb9-cd92-42ce-b724-5152216f9179" colab={"base_uri": "https://localhost:8080/", "height": 493, "referenced_widgets": ["f4b548eb2517483d9e42553b926014e5", "adcac12c6ecf4e78bb65ff2c0139b744", "b4429e47494c4ea99017f7d81c0c3320", "dadb47a0fe844961b38a07f38b047bfc", "1bdf105e4196445d9e68164039269920", "a95399e471214664b3eec7a83d67852f", "<KEY>", "<KEY>"]}
import os.path

grouped = massachusetts.groupby(['tract', 'county'])
tracts = grouped[['tract', 'county']]
tract_areas = {}
def extract_tract_area(tract):
  (this_tract, this_county) = tract[0]
  geoid = '25'+this_county+this_tract
  fc = ee.FeatureCollection("TIGER/2010/Tracts_DP1").filter(ee.Filter.eq('geoid10', geoid))
  test = fc.getInfo()['features']
  feature = ee.Feature(test[0])
  aoi = feature.geometry()
  tract_areas[geoid] = aoi.area().divide(1000 * 1000).getInfo();


for tract in tqdm(tracts):
  extract_tract_area(tract)

# + [markdown] id="O4HBOcEDAPlX" colab_type="text"
# To prevent having to run this again, let me pickle it.
#

# + id="PO4FzkGVAMge" colab_type="code" colab={}
# import pickle
# pickle.dump( tract_areas, open( "/content/drive/My Drive/areas.p", "wb" ) )

# + id="Njl0f9cQATnY" colab_type="code" colab={}
import pickle

tract_areas = favorite_color = pickle.load( open( "/content/drive/My Drive/AM216/geospatial_analysis/data/areas.p", "rb" ))


# + [markdown] id="QdhEg0XSIPwE" colab_type="text"
# Another issue is the images are not the same size, though their largest dimension is 512. For example, some may be 459x512, or 500x512. We thus create a function to pad the images like so:

# + id="lFZtG4C_AZt1" colab_type="code" colab={}
def pad_images(x):
  padded_x = np.zeros((512,512,3), dtype=np.uint8)
  x0 = x

  if x0.shape[0] != 512:
    pad_left = (512 - x0.shape[0])//2
    pad_right = x0.shape[0]+pad_left
    padded_x[pad_left:pad_right, :x0.shape[1],:] = x0

  if x0.shape[1] != 512:
    pad_top = (512 - x0.shape[1])//2
    pad_bottom = x0.shape[1]+pad_top
    padded_x[:x0.shape[0], pad_top:pad_bottom,:] = x0
  return np.array(padded_x)
  


# + [markdown] id="wFfFV5IFIdqy" colab_type="text"
# Here is a sample of the padded images:

# + id="DwLN4hJNI0pv" colab_type="code" outputId="a7750053-e181-4dfe-cd52-489c5f9ff373" colab={"base_uri": "https://localhost:8080/", "height": 286}
import imageio
from tqdm.auto import tqdm
import matplotlib.pyplot as plt
im = imageio.imread('/content/drive/My Drive/summer_resized/25001010100.png')
plt.imshow(im)

# + [markdown] id="oVSpOY1VJewp" colab_type="text"
# ### Tensorflow Data
# Given the size of the images, storing all of them in memory is not feasible. We thus leverage Tensorflow's `ImageDataGenerator` and `flow_from_dataframe` methods, which allow us to use an auxiliary Pandas dataframe containing the filename of the image as well as the population figure.
#
# First we will create this DataFrame, using CenPy

# + id="hgx_BaM6OIj-" colab_type="code" outputId="63aa5aa9-eeb0-462a-bc47-e998429b63a0" colab={"base_uri": "https://localhost:8080/", "height": 120, "referenced_widgets": ["0f3b0c766a04410f9e0d3246b88f482c", "9e9d2fd3955b4d2cb0730282d9a3ab8e", "230035818853492aabd20ebb7123c28a", "61ea164ffc064f2196f44b192e1caf42", "ef2f8d2625614022b852aba708723157", "bfd6abe99b1e41459159803c38ec0622", "8c5ab01583d4452ea6331672b38fb043", "afc364a46c2c492082f46a7b077afba4"]}
import numpy as np
from tqdm.auto import tqdm
import pandas as pd
import h5py
from matplotlib import pyplot as plt
from sklearn.model_selection import train_test_split
import pathlib

import folium
import tensorflow as tf
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import Dense, Dropout, Flatten, Activation, Input, GlobalAveragePooling2D
from tensorflow.keras.layers import Conv2D, MaxPooling2D, UpSampling2D, Concatenate
from tensorflow.keras.optimizers import Adam, SGD
from matplotlib import pyplot as plt

import os
from tensorflow.python.platform import gfile
from tensorflow.python.lib.io import file_io
from io import BytesIO
data_dir = pathlib.Path('/content/drive/My Drive/summer_resized/')

# list_ds = tf.data.Dataset.list_files(str(data_dir/'*/*'))
image_count = len(list(data_dir.glob('*.png')))
print(image_count)
IMAGE_NAMES = np.array([item.name for item in data_dir.glob('*.png') if item.name != "LICENSE.txt"])
print(IMAGE_NAMES)
POPULATION_VALS = []
for filename in tqdm(IMAGE_NAMES):
  state = filename[:2]
  county = filename[2:5]
  tract = filename[5:-4]
  pop = massachusetts[(massachusetts['state'] == state) & (massachusetts['county'] == county) & (massachusetts['tract'] == tract)].P001001.values[0]
  POPULATION_VALS.append((filename,pop))
pop_df = pd.DataFrame(POPULATION_VALS, columns =['filename', 'population']) 



# + id="kIH7x9TEZZpg" colab_type="code" colab={}
BATCH_SIZE = 16
IMG_HEIGHT = 256
IMG_WIDTH = 256
image_dir = '/content/drive/My Drive/summer_resized'

STEPS_PER_EPOCH = np.ceil(image_count/BATCH_SIZE)


# + id="gmb_RU_-xkbn" colab_type="code" outputId="d03a53f8-b221-4b66-ff20-d909e72b5042" colab={"base_uri": "https://localhost:8080/", "height": 125}
from google.colab import drive
drive.mount('/content/drive')

# + [markdown] id="dvKrruAPNmbC" colab_type="text"
# Now create the ImageDataGeneator, which will allow us to flow images.

# + id="3gt8IBe-JkDJ" colab_type="code" colab={}
train_datagen = tf.keras.preprocessing.image.ImageDataGenerator(rescale = 1./255, validation_split=0.2) 

# + [markdown] id="0iRQY02gZ-F8" colab_type="text"
# Finally, we create the generators for the training and validation.

# + id="uYk3lmPTaG2d" colab_type="code" outputId="efdbb4aa-ffec-4c3c-deab-c6e699016384" colab={"base_uri": "https://localhost:8080/", "height": 34}
train_generator = train_datagen.flow_from_dataframe(dataframe=pop_df, directory=image_dir, 
                                              x_col="filename", y_col="population", has_ext=True, 
                                              class_mode="other", target_size=(IMG_WIDTH, IMG_HEIGHT), 
                                              shuffle=True,
                                              batch_size=BATCH_SIZE,subset='training')

# + id="mLMNMj0yaOe8" colab_type="code" outputId="c78a4d75-a4ea-4b69-a353-5e0d1a3f7486" colab={"base_uri": "https://localhost:8080/", "height": 34}
 validation_generator = train_datagen.flow_from_dataframe(dataframe=pop_df, directory=image_dir, 
                                              x_col="filename", y_col="population", has_ext=True, 
                                              class_mode="other", target_size=(IMG_WIDTH, IMG_HEIGHT), 
                                              batch_size=BATCH_SIZE,subset='validation')


# + [markdown] id="mFpeFvVjvX_L" colab_type="text"
# # Models
# **Note:** The weights for our models are on our github repository, found [here](https://github.com/chickert/geospatial_analysis).
#
#

# + [markdown] id="-z9xNbue4yCk" colab_type="text"
# ## FFNN
# Our first model was a simple feed-forward neural network. This models flattens the image and then has two dense layers, before ouputting a number using linear activation. Note the summary below
#

# + id="u9bhXxn7atv8" colab_type="code" outputId="cc91e6a3-80a5-438c-84cd-3d8b17472800" colab={"base_uri": "https://localhost:8080/", "height": 336}
input1 = Input(shape=(512,512,3))
flat = Flatten()(input1)
dense1 = Dense(32, activation='relu')(flat)
dense2 = Dense(32, activation='relu')(dense1)
output = Dense(1, activation='linear')(dense2)
FFNN = Model(inputs=input1, outputs=output)
FFNN.summary()

# + id="L23CBZcRzLqY" colab_type="code" colab={}
STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
FFNN.compile(loss='mean_squared_error', optimizer='adam')


# + [markdown] id="8CI-iRSizjJu" colab_type="text"
# Now we'll train it for 10 epochs, using MSE as loss and `adam` as the optimizer.

# + id="c7IEXTvrzVJX" colab_type="code" outputId="b0339f29-2020-4d06-a255-9bc7d159a1b4" colab={"base_uri": "https://localhost:8080/", "height": 52}
# model_history = FFNN.fit(train_generator, steps_per_epoch=STEP_SIZE_TRAIN,validation_data = validation_generator, validation_steps=STEP_SIZE_VALID, epochs = 10)
FFNN.load_weights('/content/drive/My Drive/AM216/geospatial_analysis/models/FFNN.h5')
FFNN.evaluate(validation_generator, steps=STEP_SIZE_VALID)

# + [markdown] id="fXWFaNgE3x5Q" colab_type="text"
# Above we've evlauated this model on our validation dataset. The validation of loss of ~3752094 is, frankly, bad. To put this in more relatable terms, this means that each estimate of population was, on average, off by 1,937 people. We will thus attempt to improve on this with CNNs.

# + [markdown] id="XWiQTiAHBbRj" colab_type="text"
# ## Simpler CNN Models
#
# The main motivation behind buidling simpler, shallower CNN models was because we wanted to test how well models with CNN layers performed in comparison to our FFNN model with only dense layers. Generally, convolutional neural networks are able to represent image data better than dense layers because CNNs can better capture the spatial and temporal dependencies in an image.
#
# Note that due to the large amount of data, we found that adding more filters or more layers would require too much memory than google colab, could allocate (OOM error). This is another reason we started with building smaller models than our deep sequential or deep unet - we decided train them for longer (30 epochs). We see the architecture for these models below. Note: we will run the unet models on an AWS instance - details are included below.
#
# Reference source for CNNs: https://towardsdatascience.com/a-comprehensive-guide-to-convolutional-neural-networks-the-eli5-way-3bd2b1164a53

# + [markdown] id="nfMQelavB0Ob" colab_type="text"
# ### CNN Model 1
#
# - 2 Conv2D layers with filter size = 64, and relu activation
# - 2 Conv2D layers with filter size = 64, and relu activation
# - Batch normalization layers in between Conv2D layers 
# - Maxpooling in between Conv2D layers 
# - 1 Dense layer with 32 filters, relu activation
# - 1 Dropout layer
# - 1 output dense layer with linear activation
#
# Note, we experimented with many learning rates and found that 0.0001 resulted in optimal model performance. 

# + id="p7URKIAEBxf3" colab_type="code" outputId="2ebef0a1-8636-4a7f-a8de-288e5f7f2312" colab={"base_uri": "https://localhost:8080/", "height": 638}
BATCH_SIZE = 8
IMG_HEIGHT = 512
IMG_WIDTH = 512
STEPS_PER_EPOCH = np.ceil(image_count/BATCH_SIZE)
train_datagen = tf.keras.preprocessing.image.ImageDataGenerator(rescale = 1./255, validation_split=0.2)
train_generator = train_datagen.flow_from_dataframe(dataframe=pop_df, directory=image_dir,
                                              x_col="filename", y_col="population", has_ext=True,
                                              class_mode="other", target_size=(IMG_WIDTH, IMG_HEIGHT),
                                              batch_size=BATCH_SIZE,subset='training')
validation_generator = train_datagen.flow_from_dataframe(dataframe=pop_df, directory=image_dir,
                                              x_col="filename", y_col="population", has_ext=True,
                                              class_mode="other", target_size=(IMG_WIDTH, IMG_HEIGHT),
                                              batch_size=BATCH_SIZE,subset='validation')
# build model
cnn_model1 = Sequential()
cnn_model1.add(Input(shape=(512,512,3)))
cnn_model1.add(Conv2D(64,(3,3), padding="same", activation='relu'))
cnn_model1.add(tf.keras.layers.BatchNormalization())
cnn_model1.add(Conv2D(64,(3,3), padding="same", activation='relu'))
cnn_model1.add(MaxPooling2D((2,2), strides=(2,2)))
cnn_model1.add(Conv2D(32,(3,3), padding="same", activation='relu'))
cnn_model1.add(tf.keras.layers.BatchNormalization())
cnn_model1.add(Conv2D(32,(3,3),padding="same", activation='relu'))
cnn_model1.add(MaxPooling2D((2,2), strides=(2,2)))
cnn_model1.add(Flatten())
cnn_model1.add(Dense(32, activation='relu'))
cnn_model1.add(Dropout(0.3))
cnn_model1.add(Dense(1, activation='linear'))
print(cnn_model1.summary())

# build model
STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
cnn_model1.compile(loss='mean_squared_error', optimizer=Adam(0.0001))

# + id="o3Hk5WbXGdEf" colab_type="code" colab={}
# cnn_model1_history = cnn_model1.fit(train_generator,
#                           steps_per_epoch=STEP_SIZE_TRAIN,
#                           validation_data = validation_generator,
#                           validation_steps=STEP_SIZE_VALID,
#                           epochs = 30)
# cnn_model1.save_weights("model3_weights.h5")
cnn_model1.load_weights('/content/drive/My Drive/AM216/geospatial_analysis/models/cnn1.h5')
cnn_model1.evaluate(validation_generator, steps=STEP_SIZE_VALID)

# + [markdown] id="Gp_wZLJqG2Ha" colab_type="text"
# Given the simplicity of this model, it performed decently, since its best validation MSE score was 2012269.8333, which implies that the model had an average population error of 1418.54. This is better performance than the FFNN average error, although it can still be improved. We also tried implementing a simpler CNN model with CNN layers that have 32 filters.

# + [markdown] id="ie0uzHZAG6Z-" colab_type="text"
# ### CNN Model 2
#
# - 2 Conv2D layers with filter size = 32, and relu activation
# - 1 Batch normalization layer in between Conv2D layers 
# - 1 Dense layer with 32 filters, relu activation
# - 1 Dropout layer
# - 1 output dense layer with linear activation
#
# Note, we experimented with many learning rates and found that 0.00001 resulted in optimal model performance. 

# + id="nCCy5xgTG9LB" colab_type="code" outputId="b43431b8-c63b-4efe-954e-a145254bb43d" colab={"base_uri": "https://localhost:8080/", "height": 425}
cnn_model2 = Sequential()
cnn_model2.add(Input(shape=(512,512,3)))
cnn_model2.add(Conv2D(32,(3,3), padding="same", activation='relu'))
cnn_model2.add(tf.keras.layers.BatchNormalization())
cnn_model2.add(Conv2D(32,(3,3), padding="same", activation='relu'))
cnn_model2.add(Flatten())
cnn_model2.add(Dense(32, activation='relu'))
cnn_model2.add(Dropout(0.3))
cnn_model2.add(Dense(1, activation='linear'))
print(cnn_model2.summary())

# compile
STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
cnn_model2.compile(loss='mean_squared_error', optimizer=Adam(0.00001))

# + id="RPmvd66aHAEh" colab_type="code" colab={}
# # fit
# # cnn_model2_history = cnn_model2.fit(train_generator,
# #                           steps_per_epoch=STEP_SIZE_TRAIN,
# #                           validation_data = validation_generator,
# #                           validation_steps=STEP_SIZE_VALID,
# #                           epochs = 30)
# # cnn_model2.save_weights("model4_weights.h5")
cnn_model2.load_weights('/content/drive/My Drive/AM216/geospatial_analysis/models/cnn2.h5')
cnn_model2.evaluate(validation_generator, steps=STEP_SIZE_VALID)

# + [markdown] id="v0cLUQVmHG1H" colab_type="text"
# We also tried to decrease filtersize, and run for many epochs to see if this would improve performance. The lowest validation MSE loss achieved was 2531170.9358, which means the model had an average population error of 1590.96. This did not perform as well as model3, but it still is better than the FFNN. We still would like to improve performance by experimenting with different architectures and depths, so we decided to move our model fitting to AWS that allocates enough memory and speeds up performance of our models.

# + [markdown] id="Cqmm9B4BZOFi" colab_type="text"
# ## Deep Sequential CNN
#
# Given the improvement in performance we often see by deepening networks that are already partially successful, we investigated the success of deeper sequential CNNs. This investigation introduced new issues. The Google Colab notebook sessions were liable to crash due to out-of-memory (OOM) errors, and even when they did not crash, the training was exceptionally slow. 
#
# To address this, we migrated the code to AWS EC2 virtual machines, using the p3.2xlarge instance and the Anaconda environments that come pre-installed with the Deep Learning Amazon Machine Image. These instances use Tesla V100 GPUs, which come with 16GB of RAM and were sufficient to prevent OOM issues. They also trained the models much more quickly, allowing us to iterate on different design features more quickly for our comparisons. 
#
#
# This model achieved a validation loss as low as ~2.16mn, which amounts to miscounting the population of a census tract by an average of 1,469 people. Further motivation behind our specific architecture can be found [here](http://web.stanford.edu/class/archive/cs/cs221/cs221.1192/2018/restricted/posters/kfuhs/poster.pdf).
#
#
#

# + id="_h6Mpgy_ZLrb" colab_type="code" outputId="46fb0ee2-e097-4ec5-d04a-2207ba1bd9e8" colab={"base_uri": "https://localhost:8080/", "height": 816}
print("building model")
deep_model = Sequential()
deep_model.add(Input(shape=(512,512,3)))
deep_model.add(Conv2D(64,(3,3), padding="same", activation='relu'))
deep_model.add(Conv2D(64,(3,3), padding="same", activation='relu'))
deep_model.add(MaxPooling2D((2,2), strides=(2,2)))
deep_model.add(Conv2D(128,(3,3), padding="same", activation='relu'))
deep_model.add(Conv2D(128,(3,3),padding="same", activation='relu'))
deep_model.add(MaxPooling2D((2,2), strides=(2,2)))
deep_model.add(Conv2D(256,(3,3), padding="same", activation='relu'))
deep_model.add(Conv2D(256,(3,3),padding="same", activation='relu'))

deep_model.add(UpSampling2D((2,2)))
deep_model.add(Conv2D(128,(3,3), padding="same", activation='relu'))
deep_model.add(Conv2D(128,(3,3),padding="same", activation='relu'))
deep_model.add(UpSampling2D((2,2)))
deep_model.add(Conv2D(64,(3,3), padding="same", activation='relu'))
deep_model.add(Conv2D(64,(3,3),padding="same", activation='relu'))
deep_model.add(Conv2D(1,(1,1),padding="same", activation='linear'))
deep_model.add(Flatten())
deep_model.add(Dropout(0.5))
deep_model.add(Dense(1, activation='linear')) 

deep_model.summary()

STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
optimizer = Adam()
deep_model.compile(loss='mean_squared_error', optimizer=optimizer)

# fit
# model_history = deep_model.fit(train_generator, 
#                           steps_per_epoch=STEP_SIZE_TRAIN,
#                           validation_data = validation_generator, 
#                           validation_steps=STEP_SIZE_VALID, 
#                           epochs = 12)

# deep_model.save_weights("./deepseqcnn_weights.h5")

# + id="EylJEtUrZfFI" colab_type="code" colab={}
deep_model.load_weights("./deepseqcnn_weights.h5")

# + [markdown] id="7DIzC9XRabSY" colab_type="text"
# ## UNet Model
#
# To see if we could improve upon our result using the deep sequential CNN with the learning rate scheduler, we next experimented with the UNet model. This model was initially created for image segmentation problems, where it outputs a complete high resolution image. In contrast to the deep sequential model above, UNet's decoder contains skip connection layers, which combine the output of feature maps of the encoder with previous convolution layers. We can tailor
# this model to our problem, where the model can potentially classify residential areas vs non-residential areas these satellite images. Then, because our problem is a regression,
# we added layers afterwards to output a parameter representing population count. To reduce problem size and overfitting, we also added a dropout layer in between the two dense layers. More information on the [original UNet paper](https://arxiv.org/pdf/1505.04597.pdf) for reference, and further motivation behind building this model can be found [here](http://web.stanford.edu/class/archive/cs/cs221/cs221.1192/2018/restricted/posters/kfuhs/poster.pdf).
#
# This resulted in a validation loss as low as ~2.28mn, which amounts to miscounting population by an average of 1,509 per census tract. 
#
#

# + id="i9XkC2eHaeZD" colab_type="code" outputId="0bbe4c5a-6b5e-41c4-f2d3-66ee1559e0c0" colab={"base_uri": "https://localhost:8080/", "height": 976}
print("building model")
input_feat = Input(shape=(512,512,3))
conv_layer1 = Conv2D(64,(3,3), padding="same", activation='relu')(input_feat)
conv_layer2 = Conv2D(64,(3,3), padding="same", activation='relu')(conv_layer1)
max_pool1 = MaxPooling2D((2,2), strides=(2,2))(conv_layer1)
conv_layer3 = Conv2D(128,(3,3), padding="same", activation='relu')(max_pool1)
conv_layer4 = Conv2D(128,(3,3), padding="same", activation='relu')(conv_layer3)
max_pool2 = MaxPooling2D((2,2), strides=(2,2))(conv_layer4)
conv_layer5 = Conv2D(256,(3,3), padding="same", activation='relu')(max_pool2)
conv_layer6 = Conv2D(256,(3,3),padding="same", activation='relu')(conv_layer5)

# upsampling and concatenating
upsample1 = UpSampling2D((2,2))(conv_layer6)
concat1 = Concatenate()([upsample1, conv_layer4])
conv_layer7 = Conv2D(128,(3,3), padding="same", activation='relu')(concat1)
conv_layer8 = Conv2D(128,(3,3),padding="same", activation='relu')(conv_layer7)
upsample2 = UpSampling2D((2,2))(conv_layer8)
concat2 = Concatenate()([upsample2, conv_layer2])
conv_layer9 = Conv2D(64,(3,3), padding="same", activation='relu')(concat2)
conv_layer10 = Conv2D(64,(3,3),padding="same", activation='relu')(conv_layer9)
output = Conv2D(1,(1,1),padding="same", activation='linear')(conv_layer10)

reshape_layer1 = Flatten()(output)
reshape_layer2 = Dropout(0.3)(reshape_layer1)
final_output = Dense(1, activation='linear')(reshape_layer2)

# build model
unet_model = tf.keras.Model(input_feat, final_output)
unet_model.summary()

STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
unet_model.compile(loss='mean_squared_error', optimizer='adam')

# fit
print("training...")
# model_history = unet_model.fit(train_generator,
#                           steps_per_epoch=STEP_SIZE_TRAIN,
#                           validation_data = validation_generator,
#                           validation_steps=STEP_SIZE_VALID,
#                           epochs = 12)

# unet_model.save_weights("./unet_weights.h5")
unet_model.load_weights("./unet_weights.h5")

# + [markdown] id="e8jULD1Vamym" colab_type="text"
# ## UNet Model with Dense Layers
#
# The UNet archictecture was originally designed for image segmentation with the ultimate purpose of classification. However, since our use-case is a regression problem, we experimented with adding several dense layers (in addition to the one in the model above necessary to output a population estimate) after the convolutional segment of the network.
#
# However, we found that this did not improve performance. The model architecture is defined below, and only drops the validation loss as low as ~2.62mn, which amounts to miscounting the population by 1,618 heads per census tract. 
#

# + id="X-aq8MNDaris" colab_type="code" outputId="3ca9a71b-7e97-49ca-e10f-a0f8158a7a94" colab={"base_uri": "https://localhost:8080/", "height": 1000}
print("building model")
input_feat = Input(shape=(512,512,3))
conv_layer1 = Conv2D(64,(3,3), padding="same", activation='relu')(input_feat)
conv_layer2 = Conv2D(64,(3,3), padding="same", activation='relu')(conv_layer1)
max_pool1 = MaxPooling2D((2,2), strides=(2,2))(conv_layer1)
conv_layer3 = Conv2D(128,(3,3), padding="same", activation='relu')(max_pool1)
conv_layer4 = Conv2D(128,(3,3), padding="same", activation='relu')(conv_layer3)
max_pool2 = MaxPooling2D((2,2), strides=(2,2))(conv_layer4)
conv_layer5 = Conv2D(256,(3,3), padding="same", activation='relu')(max_pool2)
conv_layer6 = Conv2D(256,(3,3),padding="same", activation='relu')(conv_layer5)

# upsampling and concatenating
upsample1 = UpSampling2D((2,2))(conv_layer6)
concat1 = Concatenate()([upsample1, conv_layer4])
conv_layer7 = Conv2D(128,(3,3), padding="same", activation='relu')(concat1)
conv_layer8 = Conv2D(128,(3,3),padding="same", activation='relu')(conv_layer7)
upsample2 = UpSampling2D((2,2))(conv_layer8)
concat2 = Concatenate()([upsample2, conv_layer2])
conv_layer9 = Conv2D(64,(3,3), padding="same", activation='relu')(concat2)
conv_layer10 = Conv2D(64,(3,3),padding="same", activation='relu')(conv_layer9)
output = Conv2D(1,(1,1),padding="same", activation='linear')(conv_layer10)

added_layer1 = Flatten()(output)
added_layer2 = Dropout(0.3)(added_layer1)
added_layer3 = Dense(256, activation='relu')(added_layer2)
added_layer4 = Dropout(0.3)(added_layer3)
added_layer5 = Dense(64, activation='relu')(added_layer4)
added_layer6 = Dropout(0.3)(added_layer5)
final_output = Dense(1, activation='linear')(added_layer6)

# build model
unet_model = tf.keras.Model(input_feat, final_output)
unet_model.summary()

STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
unet_model.compile(loss='mean_squared_error', optimizer='adam')

# fit
# print("training...")
# model_history = unet_model.fit(train_generator,
#                           steps_per_epoch=STEP_SIZE_TRAIN,
#                           validation_data = validation_generator,
#                           validation_steps=STEP_SIZE_VALID,
#                           epochs = 11)

# unet_model.save_weights("./unetwdenselayers_weights.h5")
unet_model.load_weights("./unetwdenselayers_weights.h5")

# + [markdown] id="XO6RkalU57eb" colab_type="text"
# ## Transfer Learning
# While the task of population prediction is relatively distinct from any sort of classification task, we created a model build on top of MobileNet with a few dense layers added to the end before the linear activation layer. More information about MobileNet is availible here: https://arxiv.org/abs/1704.04861. Considering this was a last-ditch effort, this model was a pleasant surprise. Below I define the model and evaluate it on our validation dataset.

# + id="zGpl6_4VzlnG" colab_type="code" outputId="484122c3-94b9-4a69-d936-ed15e164f701" colab={"base_uri": "https://localhost:8080/", "height": 1000}
from tensorflow.keras.applications import MobileNet

base_model=MobileNet(weights='imagenet',include_top=False, input_shape=(256,256,3)) #imports the mobilenet model and discards the last 1000 neuron layer.

x=base_model.output
x=GlobalAveragePooling2D()(x)
x=Dense(1024,activation='relu')(x) #we add dense layers so that the model can learn more complex functions and classify for better results.
x=Dense(1024,activation='relu')(x) #dense layer 2
x=Dense(512,activation='relu')(x) #dense layer 3
preds=Dense(1,activation='linear')(x) 
transfer_model = Model(base_model.input, preds)
transfer_model.summary()

# + [markdown] id="5TSPO5HvD4Kh" colab_type="text"
# Now we evaluate it on the validation

# + id="1gcqd5ZsD_HH" colab_type="code" outputId="790c7324-c23a-4ed5-ad2e-42eb301ca9fd" colab={"base_uri": "https://localhost:8080/", "height": 260}
STEP_SIZE_TRAIN=train_generator.n//train_generator.batch_size
STEP_SIZE_VALID=validation_generator.n//validation_generator.batch_size
transfer_model.compile(loss='mean_squared_error', optimizer='adam')
transfer_model.load_weights('/content/drive/My Drive/AM216/geospatial_analysis/models/trans.h5')
# model_history = transfer_model.fit(train_generator, steps_per_epoch=STEP_SIZE_TRAIN,validation_data = validation_generator, validation_steps=STEP_SIZE_VALID, epochs = 10)

# transfer_model.evaluate(validation_generator, steps=STEP_SIZE_VALID)

# + [markdown] id="nOqugRSC14w1" colab_type="text"
# As you can see above, the validation loss of ~1871448 equates an average error in population of 1368. While certainly not perfect, this is a marked improvement over many of our other models. This is a testement to the power of transfer learning.

# + [markdown] id="F-o_5D3CuPk8" colab_type="text"
# # Diving into our best model
# Let's quickly take a look at the outcomes of our predictions, using the Transfer Learning model.
#
# Here's the loss history plot for our model. It seems like a high degree of overfitting is going on.
#

# + id="FRQIeEDK0UI9" colab_type="code" outputId="81e2a156-161f-43d5-c7dd-b233e580f5be" colab={"base_uri": "https://localhost:8080/", "height": 286}
# import imageio
from tqdm.auto import tqdm
import matplotlib.pyplot as plt
im = imageio.imread('/content/drive/My Drive/AM216/geospatial_analysis/history.png')
plt.imshow(im)

# + id="aco4nvPEujSd" colab_type="code" outputId="0390dcdd-96e8-4c42-99c0-9df5bde181ff" colab={"base_uri": "https://localhost:8080/", "height": 66, "referenced_widgets": ["78cad7b29fe8413e9f61fe8de0d7b39a", "9124e0d7aaf4416ea5a6feeb254ba6f0", "28fcc93993f548c5910b66a1cc0d058c", "81ead2f87d8247f49a9b243256654e0d", "e8554fd887bf40c1b721494a1ed55c7f", "cd9bc1b656124c89b021ce139a2110bd", "2fb40b80285a4285b3339a484c8f913e", "8c444a4291b843f4940963b2d26274b5"]}
max_diff = 0;
real_pop = 0
predicted_pop = 0 
for i in tqdm(range(18)):
    x,y = next(train_generator)
    preds = transfer_model.predict(x)
    diff = np.abs(preds.flatten() - y)
    local_max = np.max(diff)
    if local_max > max_diff:
        max_diff = local_max
        real_pop = y[np.argmax(diff)]
        predicted_pop = preds[np.argmax(diff)]

# + id="maEvBV1V1hCk" colab_type="code" outputId="abbc1db5-23a0-45a4-d804-fca4b7a61085" colab={"base_uri": "https://localhost:8080/", "height": 77}
pop_df[pop_df.population == 6652.0]


# + id="yMQobmj31zeQ" colab_type="code" outputId="2f5c2b06-414d-4d47-ebbb-ffbeb3982052" colab={"base_uri": "https://localhost:8080/", "height": 286}
im = imageio.imread('/content/drive/My Drive/summer_resized/25017340100.png')
plt.imshow(im)

# + [markdown] id="MRpQZQkR_BNx" colab_type="text"
# Here's the tract that our model performed the worst on. It looks like it's a large, suburban neighborhood. Perhaps the large amount of green space confused the model, and if we had a larger dataset and/or higher definition photos, the model would behave better. 

# + [markdown] id="dc5M-yPW5BHD" colab_type="text"
# # Unsuccessful Forays
# We made a number of attempts to improve the performance of our models.

# + [markdown] id="gfFXM8dD5UAF" colab_type="text"
# ## Adding area of tract as a second input
# We thought that perhaps an issue with our data was the fact that each image of a tract was the same size, despite tracts ranging widely in size. Thus we created a new data generator which combined the image as well as the area of the tract. Below is the code for creating such a generator.
#
#

# + id="CvW3C1qG5ONY" colab_type="code" colab={}
train_datagen = tf.keras.preprocessing.image.ImageDataGenerator(rescale = 1./255, validation_split=0.2) 
def combined_generator(subset, width, height, batchsize):
    img_gen = train_datagen.flow_from_dataframe(dataframe=pop_df, directory=image_dir, 
                                              x_col="filename", y_col="population", has_ext=True, 
                                              class_mode="other", target_size=(width, height), 
                                              batch_size=batchsize,subset=subset, shuffle=True)
    while True:
      idx = (img_gen.batch_index) * img_gen.batch_size
      filenames = img_gen.filenames[idx : idx + img_gen.batch_size]
      area_list = pop_df[pop_df['filename'].isin(filenames)].area.values
      X1 = next(img_gen)
      yield ([X1[0],area_list], X1[1])
train_generator = combined_generator(subset='training', width=IMG_WIDTH, height=IMG_HEIGHT, batchsize=BATCH_SIZE)
validation_generator = combined_generator(subset='validation', width=IMG_WIDTH, height=IMG_HEIGHT, batchsize=BATCH_SIZE)


# + [markdown] id="F9mnrx5lpzoH" colab_type="text"
# ## Adding a learning rate scheduler to fit the models
#
# We experimented with adding a learning rate scheduler to improver model performance. A learning scheduler audo-adjusts the learninng rate per epoch to performance at each epoch of training. For us, the learning rate scheduler started with a train learning rate of 0.001, and validation learning rate of 0.0001. we found that although it improved train MSE, it did not improve validation MSE at the same rate, resulting in overfitting. As a result, we left it out of our final models. You can incorporate a learning rate scheduler as demonstrated below.

# + id="lC-3e-Kjp0Sx" colab_type="code" colab={}
# functions for learning rate scheduler
def get_lr_metric(optimizer): # to print out learning rate during model fitting
    def lr(y_true, y_pred):
        return optimizer.lr
    return lr

def scheduler(epoch): # exponential decay of learning rate
  initial_lrate = 0.0001
  if epoch < 5:
    lrate = initial_lrate
  else:
    lrate = initial_lrate*math.exp(0.1*(5-epoch))
  return lrate

# compile desired model
optimizer = Adam()
lr_metric = get_lr_metric(optimizer)
unet_model.compile(loss='mean_squared_error', optimizer=optimizer, metrics=['accuracy', lr_metric])

# Define learning rate scheduler
lrate = LearningRateScheduler(scheduler)
callbacks_list = [lrate]

# Fit
unet_model_history = unet_model.fit(train_generator,
                          steps_per_epoch=STEP_SIZE_TRAIN,
                          validation_data = validation_generator,
                          validation_steps=STEP_SIZE_VALID, callbacks=callbacks_list, epochs=15)

# + [markdown] id="TT8SHIl05yj0" colab_type="text"
# # Discussion of Results
# ### Potential Pitfalls
# A number of reasons may have hindered us. The highest resolution satellite image data is not free, and has more than 10x the resolution of the data that we used. In using this DigitalGlobe data, we could have conducted our analyses at the smaller "census block" level. Perhaps this would have been more successful, and the models would have learned more abou the nature of human habitation.
#
# Given the time it took to train our models, we limited ourselves to data from Massachusetts, as opposed to the entire United States. It would have been more interesting to see how our approach fared on a broader dataset, and perhaps our models would have done better in this case, as there is only so much a model can learn from fewer than 2000 total images. 
#
# Hence, computational power was a significant limitation, and though we attempted to get around it by spinning up GPU EC2 instances as opposed to training within colab, we still did not have access large, multi-GPU instance, which could have allowed us to look at more data.
#
# Finally, there is the foundational difficulties that come with the problem of predicitng populations based on landscape images alone. For example, intuitively we might think of buildings as correlating with higher population density. This may be true in many places such as suburban neighborhoods or areas with dense apartment complex clusters, but financial districts in city centers may -- on the contrary -- be replete with buildings, but have near-zero population counts in the census. Put another way, what features readily distinguish an apartment skyscraper from a commercial office space skyscraper? 
#
# ### FFNN
# We created a FFNN with only a few dense layers as a basleine model to compare our other prediction results to. We found that the model had a poor validation of loss of ~3752094, or in other words, each estimate of population was, on average, off by 1,937 people. We anticipated this model to perform poorly, as it was not tweaked for optimal performance.
#
# ### Simple CNN models
#
# We found that the CNN models with fewer layers were able to predict population with less average population count error per census track than FFNN (average error of ~1500 people). We were surprised that the average error wasn't too far off, given the simple CNN archtiecture, and that the task of detecting highly populated areas purely from satellite images is relatively difficult. There were a decent amount of satellite images without much population (i.e. Western MA), and so we think that the model would more easily tell the difference between populated vs unpopulated areas. We think that the CNN models with less layers has a limit with how accuractely they could predict (could not achieve average population area less than ~1500 people), and so these models were limited. 
#
# ### Deep Sequential Model and UNet Models: 
#
# The deep sequential and UNet models both consisted of more layers. The difference between the two is that the UNet model's decoder contains skip connection layers, which combine the output of feature maps of the encoder with previous convolution layers. We were disappointed to see that these deeper models did not boost the performance of the simpler models, again stalling around the same ~1500 limit as described for those more shallow models.
#
# ### Transfer Learning Model
#
# The transfer learning model, built on top of MobileNet was the most performant of our models. It should still be noted that this peak performance is still an error range of +/- ~1000 people per tract, which still isn't great. However, this demonstrates the wide applicability of transfer learning, even when not attempting to do classification.
#
#
# # Conclusion
# Overall through this project, we were able to explore the GEE database, build a data pipeline to obtain the satellite images, perform data preprocssing to pad the satellite images and combine demographic data (from CenPy), and build multiple neural network models to predict population. To conclude, we were successful in building a pipeline and models that analyzed and predicted based off of satellite imagery. However, we found that training on ~1400 images was not enough for the model to learn weights that could accurately predict population. We also found that the quality of images affects performance, and that a more complex architecture (UNet model) does not necessarily outperform a simpler architecture (simple CNN models). Regardless of our models not performing incredibly well, we were able to experiment with different models and learn how to interact with GEE. We found the inverse problem of learning population from satellite images to be an interesting, challenging, and rewarding task to solve.
#

# + id="1_CGW9tc5vQS" colab_type="code" colab={}


StarcoderdataJupyterScriptsDedupFiltered