diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..3fd6e7e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*.pyc
+*.mat
\ No newline at end of file
diff --git a/box_head.py b/box_head.py
new file mode 100644
index 0000000..e6221f3
--- /dev/null
+++ b/box_head.py
@@ -0,0 +1,80 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 14:44:57 2018
+
+@author: z5030440
+
+Main code to extract shorelines from Landsat imagery
+"""
+# Preamble
+
+import ee
+import math
+import matplotlib.pyplot as plt
+import numpy as np
+import pdb
+
+# image processing modules
+import skimage.filters as filters 
+import skimage.exposure as exposure
+import skimage.transform as transform
+import sklearn.decomposition as decomposition
+import skimage.morphology as morphology
+import skimage.measure as measure
+
+# my modules
+# my functions
+from functions.utils import *
+from functions.sds import *
+
+np.seterr(all='ignore') # raise/ignore divisions by 0 and nans
+ee.Initialize()
+
+# parameters
+plot_bool = True # if you want the plots
+prob_high = 99.9 # upper probability to clip and rescale pixel intensity
+min_contour_points = 100 # minimum number of points contained in each water line
+
+# select collection
+input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
+
+# location (Narrabeen-Collaroy beach)
+rect_narra = [[[151.317395,-33.494601],
+              [151.388635,-33.495174], 
+              [151.363624,-33.565184],
+              [151.305228,-33.563299],
+              [151.317395,-33.494601]]];
+
+# filter by location
+flt_col = input_col.filterBounds(ee.Geometry.Polygon(rect_narra))
+
+n_img = flt_col.size().getInfo()
+print('Number of images covering Narrabeen:', n_img)
+im_all = flt_col.getInfo().get('features')
+output = []
+# loop through all images
+# find each image in ee database
+i = 2
+im = ee.Image(im_all[i].get('id')) 
+# load image as np.array
+im_pan, im_ms, im_cloud, crs = read_eeimage(im, rect_narra, plot_bool)
+# rescale intensities
+im_ms = rescale_image_intensity(im_ms, im_cloud, prob_high, plot_bool)
+im_pan = rescale_image_intensity(im_pan, im_cloud, prob_high, plot_bool)
+# pansharpen rgb image
+im_ms_ps = pansharpen(im_ms[:,:,[0,1,2]], im_pan, im_cloud, plot_bool)
+# add down-sized bands for NIR and SWIR (since pansharpening is not possible)
+im_ms_ps = np.append(im_ms_ps, im_ms[:,:,[3,4]], axis=2) 
+# calculate NDWI
+im_ndwi = nd_index(im_ms_ps[:,:,3], im_ms_ps[:,:,1], im_cloud, plot_bool)
+# edge detection
+wl_pix = find_wl_contours(im_ndwi, im_cloud, min_contour_points, True)
+# convert from pixels to world coordinates
+wl_coords = convert_pix2world(wl_pix, crs['crs_15m'])
+output.append(wl_coords)
+
+plt.figure()
+plt.imshow(im_ms_ps[:,:,[2,1,0]])
+plt.axis('off')
+plt.title('RGB at 15m')
+plt.show()