mirror of
https://github.com/wassname/satellite_leak_detection.git
synced 2026-09-10 12:34:13 +08:00
Draft
This commit is contained in:
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,50 @@
|
||||
|
||||
Data received from Austin Water as part of public information request: PIR 32348 on February 16, 2017. Locations are the nearest street address and no other contextual information was given.
|
||||
|
||||
Parsing the data, working out the columns and other details are in `notebooks/0_process_data/process_austin_leaks.ipynb`.
|
||||
|
||||
----------
|
||||
|
||||
Dates are in US/Central timezone.
|
||||
|
||||
Locations seem to be EPSG:6578.
|
||||
|
||||
Columns:
|
||||
|
||||
The columns are unlabelled but the top left 47486 match leaks pending inspections https://services1.arcgis.com/PuB3FWUAxkScvfQy/ArcGIS/rest/services/LPI/FeatureServer/0
|
||||
|
||||
The bottom right data matches leak pending repair https://services1.arcgis.com/PuB3FWUAxkScvfQy/ArcGIS/rest/services/LWOPR/FeatureServer/0
|
||||
|
||||
```py
|
||||
columns = [
|
||||
'OBJECTID',
|
||||
'INSPFLAG',
|
||||
'PRI',
|
||||
'PROB',
|
||||
'PROBDTTM', # Problem datetime US/Central timezone.
|
||||
'SCHEDDTTM', # Scheduled datetime US/Central timezone.
|
||||
'SERVNO',
|
||||
'MAPNO',
|
||||
'STARTDTTM', # Start datetime US/Central timezone.
|
||||
'PROBCODE',
|
||||
'PROBDESC',
|
||||
'CITY',
|
||||
'PREDIR',
|
||||
'STNAME',
|
||||
'STNO',
|
||||
'STSUB',
|
||||
'SUFFIX',
|
||||
'ZIP',
|
||||
'ADDRKEY',
|
||||
'FullStreetName',
|
||||
'X', # X location in EPSG:6578
|
||||
'Y', # Y location in EPSG:6578
|
||||
'22', # CMPLKEY (Complementary key?, INT)?, SERVER (nah thats 47486-627096, Objectid (no too high)?? 1913-1496485.
|
||||
'QTYCALLS', # Quantity of calls?
|
||||
'INITDTTM', # Initial date time US/Central timezone.
|
||||
'WONO', # Work order number?
|
||||
'LOC',
|
||||
'DESCRIPT', # Description
|
||||
'COMPDTTM' # Completion data time US/Central timezone.
|
||||
]
|
||||
```
|
||||
Binary file not shown.
@@ -0,0 +1,5 @@
|
||||
The information here is from FOI request to South Australia Water. FOI id "SN840 - FOI Determination - Leak repair work orders"
|
||||
|
||||
Date: Thursday, March 23, 2017 1:47 PM
|
||||
|
||||
The time column is report data, so hopefully it's the time of repair as I asked not the time that the leak was reported. I will clarify this...
|
||||
Binary file not shown.
File diff suppressed because one or more lines are too long
@@ -0,0 +1,16 @@
|
||||
Files:
|
||||
- script_metadata.json - information on scraping script
|
||||
- data.h5 contains X and y
|
||||
- X: tiff files for each band loaded into an array of shape (Leak, Bands, width, length)
|
||||
- y: True for before the leak, False for after
|
||||
- data_metadata: array of metadata for each leak in X. Each contain info on leak, image, and image search
|
||||
|
||||
Loading:
|
||||
```py
|
||||
# load
|
||||
metadatas = json.load(open('data_metadata.json'))
|
||||
with h5py.File('data.h5','r') as h5f:
|
||||
X2 = h5f['X'][:]
|
||||
y2 = h5f['y'][:]
|
||||
y
|
||||
```
|
||||
@@ -0,0 +1 @@
|
||||
{"temp_dir": "/tmp/scraping_earth_engine_NAIP_all_v2", "bands": ["R", "G", "B", "N"], "resolution_min": 1.0, "ts": "20170321-02-54-41", "crs_grid": 3857, "notebook_name": "scraping_earth_engine_NAIP_all_v2", "cache_dir": "../data/downloaded_images-scraping_earth_engine_NAIP_all_v2-USDA-NAIP-DOQQ/cache", "pixel_length": 129.0, "output_dir": "../data/downloaded_images-scraping_earth_engine_NAIP_all_v2-USDA-NAIP-DOQQ", "satellite": "USDA/NAIP/DOQQ", "time_bin_delta": 2419200}
|
||||
Binary file not shown.
File diff suppressed because one or more lines are too long
@@ -0,0 +1,16 @@
|
||||
Files:
|
||||
- script_metadata.json - information on scraping script
|
||||
- data.h5 contains X and y
|
||||
- X: tiff files for each band loaded into an array of shape (Leak, Bands, width, length)
|
||||
- y: True for before the leak, False for after
|
||||
- data_metadata: array of metadata for each leak in X. Each contain info on leak, image, and image search
|
||||
|
||||
Loading:
|
||||
```py
|
||||
# load
|
||||
metadatas = json.load(open('data_metadata.json'))
|
||||
with h5py.File('data.h5','r') as h5f:
|
||||
X2 = h5f['X'][:]
|
||||
y2 = h5f['y'][:]
|
||||
y
|
||||
```
|
||||
@@ -0,0 +1 @@
|
||||
{"resolution_min": 15.0, "crs_grid": 3857, "output_dir": "../data/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2", "ts": "20170318-10-08-12", "satellite": "LANDSAT/LE7_L1T", "notebook_name": "testing_earth_engine-l7-AUTX", "bands": ["B1", "B2", "B3", "B4", "B5", "B6_VCID_1", "B6_VCID_2", "B7", "B8"], "cache_dir": "../data/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2/ee_l7_AUTX-leaks_cache_v2", "pixel_length": 25.0, "temp_dir": "/tmp/testing_earth_engine-l7-AUTX-48wtrw24-20170318-10-08-12"}
|
||||
Binary file not shown.
File diff suppressed because one or more lines are too long
@@ -0,0 +1,16 @@
|
||||
Files:
|
||||
- script_metadata.json - information on scraping script
|
||||
- data.h5 contains X and y
|
||||
- X: tiff files for each band loaded into an array of shape (Leak, Bands, width, length)
|
||||
- y: True for before the leak, False for after
|
||||
- data_metadata: array of metadata for each leak in X. Each contain info on leak, image, and image search
|
||||
|
||||
Loading:
|
||||
```py
|
||||
# load
|
||||
metadatas = json.load(open('data_metadata.json'))
|
||||
with h5py.File('data.h5','r') as h5f:
|
||||
X2 = h5f['X'][:]
|
||||
y2 = h5f['y'][:]
|
||||
y
|
||||
```
|
||||
@@ -0,0 +1,13 @@
|
||||
{
|
||||
"resolution_min": 15.0,
|
||||
"notebook_name": "testing_earth_engine-l8-AUTX_v2",
|
||||
"satellite": "LANDSAT/LC8_L1T",
|
||||
"pixel_length": 25.0,
|
||||
"crs_grid": 3857,
|
||||
"time_bin_delta": 2419200,
|
||||
"cache_dir": "../data/20170314-05-26-52_testing_earth_engine-l8-AUTX_v2/ee_l8_AUTX-leaks_cache_v2/ATX_",
|
||||
"ts": "20170321-02-12-16",
|
||||
"temp_dir": "/tmp/testing_earth_engine-l8-AUTX_v2-82dzkbkv-20170321-02-12-16",
|
||||
"output_dir": "../data/20170314-05-26-52_testing_earth_engine-l8-AUTX_v2",
|
||||
"bands": ["B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B9", "B10", "B11", "BQA"]
|
||||
}
|
||||
Binary file not shown.
File diff suppressed because one or more lines are too long
@@ -0,0 +1,16 @@
|
||||
Files:
|
||||
- script_metadata.json - information on scraping script
|
||||
- data.h5 contains X and y
|
||||
- X: tiff files for each band loaded into an array of shape (Leak, Bands, width, length)
|
||||
- y: True for before the leak, False for after
|
||||
- data_metadata: array of metadata for each leak in X. Each contain info on leak, image, and image search
|
||||
|
||||
Loading:
|
||||
```py
|
||||
# load
|
||||
metadatas = json.load(open('data_metadata.json'))
|
||||
with h5py.File('data.h5','r') as h5f:
|
||||
X2 = h5f['X'][:]
|
||||
y2 = h5f['y'][:]
|
||||
y
|
||||
```
|
||||
@@ -0,0 +1 @@
|
||||
{"ts": "20170318-10-39-35", "satellite": "COPERNICUS/S1_GRD", "crs_grid": 3857, "bands": ["VV", "HH", "VH", "HV", "angle"], "resolution_min": 10.0, "output_dir": "../data/downloaded_images_scraping_earth_engine_s1-all_COPERNICUS-S1_GRD", "temp_dir": "/tmp/scraping_earth_engine_s1-all", "notebook_name": "scraping_earth_engine_s1-all", "pixel_length": 25.0, "cache_dir": "../data/downloaded_images_scraping_earth_engine_s1-all_COPERNICUS-S1_GRD/cache", "time_bin_delta": 2419200}
|
||||
Binary file not shown.
File diff suppressed because one or more lines are too long
@@ -0,0 +1,16 @@
|
||||
Files:
|
||||
- script_metadata.json - information on scraping script
|
||||
- data.h5 contains X and y
|
||||
- X: tiff files for each band loaded into an array of shape (Leak, Bands, width, length)
|
||||
- y: True for before the leak, False for after
|
||||
- data_metadata: array of metadata for each leak in X. Each contain info on leak, image, and image search
|
||||
|
||||
Loading:
|
||||
```py
|
||||
# load
|
||||
metadatas = json.load(open('data_metadata.json'))
|
||||
with h5py.File('data.h5','r') as h5f:
|
||||
X2 = h5f['X'][:]
|
||||
y2 = h5f['y'][:]
|
||||
y
|
||||
```
|
||||
@@ -0,0 +1,12 @@
|
||||
{
|
||||
"bands": ["B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B9", "B10", "B11", "B12", "QA60"],
|
||||
"crs_grid": 3857,
|
||||
"ts": "20170316-06-09-45",
|
||||
"output_dir": "../data/downloaded_images_image_testing_earth_engine_s2-AUTX_v6_COPERNICUS-S2",
|
||||
"temp_dir": "/tmp/image_testing_earth_engine_s2-AUTX_v6",
|
||||
"satellite": "COPERNICUS/S2",
|
||||
"pixel_length": 25.0,
|
||||
"resolution_min": 10.0,
|
||||
"notebook_name": "image_testing_earth_engine_s2-AUTX_v6",
|
||||
"cache_dir": "../data/downloaded_images_image_testing_earth_engine_s2-AUTX_v6_COPERNICUS-S2/cache"
|
||||
}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 59 KiB |
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -0,0 +1,108 @@
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def calculate_result_class(y_pred, y_true, thresh=0.5, words=False):
|
||||
"""Calculate int results classes. Here 'false positive' = 'fp' = bin('01') = 1 etc."""
|
||||
dic = {0: 'fn', 1: 'fp', 2: 'tn', 3: 'tp'}
|
||||
y_pred = np.array(y_pred) > thresh
|
||||
y_true = np.array(y_true) > thresh
|
||||
|
||||
result_class_bools = np.array(list(zip(y_pred == y_true, y_pred))) * 1
|
||||
# convert to binary
|
||||
result_class_binary = ['0b' + ''.join(d) for d in result_class_bools.astype(str)]
|
||||
result_class = [int(d, 2) for d in result_class_binary]
|
||||
|
||||
if words:
|
||||
result_class = [dic[d] for d in result_class]
|
||||
return result_class
|
||||
|
||||
|
||||
assert calculate_result_class([0, 1, 0, 1], [1, 0, 0, 1], words=False) == [0, 1, 2, 3]
|
||||
assert calculate_result_class([0, 1, 0, 1], [1, 0, 0, 1], words=True) == ['fn', 'fp', 'tn', 'tp']
|
||||
assert calculate_result_class([False, True, False, True], [True, False, False, True], words=False) == [0, 1, 2, 3]
|
||||
|
||||
|
||||
from io import StringIO
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
from sklearn import metrics
|
||||
|
||||
|
||||
def parse_classification_report(classification_report):
|
||||
"""Parse a sklearn classification report to a dict."""
|
||||
return pd.read_fwf(
|
||||
StringIO(classification_report),
|
||||
index_col=0,
|
||||
colspecs=[(0, 12), (12, 22), (22, 32), (32, 42), (42, 52)]
|
||||
).dropna()
|
||||
|
||||
# test
|
||||
s = metrics.classification_report(np.random.random(100) > 0.5, np.random.random(100) > 0.5)
|
||||
d = parse_classification_report(s).to_dict()
|
||||
assert isinstance(d, dict)
|
||||
|
||||
|
||||
import sklearn
|
||||
from sklearn.dummy import DummyClassifier, DummyRegressor
|
||||
from sklearn.model_selection import train_test_split
|
||||
import collections
|
||||
|
||||
|
||||
def find_best_dummy_classification(X, y, test_size=0.3, random_state=0, thresh=0.5, target_names=None, n=1):
|
||||
"""Try all dummy models."""
|
||||
X = X.reshape((len(X) ,-1))
|
||||
# y = y.reshape((len(y) ,-1))
|
||||
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=random_state)
|
||||
|
||||
dummy_scores = []
|
||||
for i in range(n):
|
||||
for strategy in ['most_frequent', 'uniform', 'prior', 'stratified']:
|
||||
clf = DummyClassifier(strategy=strategy)
|
||||
clf.fit(X_train, y_train)
|
||||
y_pred = clf.predict(X_test)
|
||||
score = clf.score(X_test, y_test)
|
||||
|
||||
matthews_corrcoef=sklearn.metrics.matthews_corrcoef(y_test > thresh, y_pred > thresh)
|
||||
|
||||
report=parse_classification_report(sklearn.metrics.classification_report(y_test > thresh, y_pred > thresh, target_names=target_names))
|
||||
|
||||
dummy_scores.append(
|
||||
collections.OrderedDict(
|
||||
strategy='classifier_' + strategy,
|
||||
matthews_corrcoef=matthews_corrcoef,
|
||||
score=score,
|
||||
report=report
|
||||
)
|
||||
)
|
||||
|
||||
for strategy in ['mean', 'median']:
|
||||
clf=DummyRegressor(strategy=strategy)
|
||||
clf.fit(X_train, y_train)
|
||||
y_pred=clf.predict(X_test)
|
||||
score=clf.score(X_test, y_test)
|
||||
|
||||
matthews_corrcoef=sklearn.metrics.matthews_corrcoef(y_test > thresh, y_pred > thresh)
|
||||
|
||||
report=parse_classification_report(sklearn.metrics.classification_report(y_test > thresh, y_pred > thresh, target_names=target_names))
|
||||
|
||||
dummy_scores.append(
|
||||
collections.OrderedDict(
|
||||
strategy='regressor_' + strategy,
|
||||
matthews_corrcoef=matthews_corrcoef,
|
||||
score=score,
|
||||
report=report
|
||||
)
|
||||
)
|
||||
|
||||
df=pd.DataFrame(dummy_scores)
|
||||
df=df.sort_values('matthews_corrcoef', ascending=False)
|
||||
return df, df[:1].iloc[0].to_dict()
|
||||
|
||||
# test
|
||||
# import numpy as np
|
||||
# X=np.random.random((100,10,10,10))
|
||||
# y=np.random.random((100))>0.5
|
||||
# df, dummy = find_best_dummy_classification(X,y,n=100)
|
||||
# df_dummies.groupby('strategy').mean()
|
||||
@@ -0,0 +1,254 @@
|
||||
import logging
|
||||
import numpy as np
|
||||
import ee
|
||||
import zipfile
|
||||
import tempfile
|
||||
from tqdm import tqdm
|
||||
from scipy.misc import imread
|
||||
import urllib
|
||||
from path import Path
|
||||
|
||||
ee.Initialize() # should give no errors, if so follow instructions
|
||||
logger = logging.getLogger('leaks_helpers')
|
||||
temp_dir = Path(tempfile.mkdtemp())
|
||||
|
||||
bands_NAIP = ['R', 'G', 'B', 'N']
|
||||
|
||||
# params https://explorer.earthengine.google.com/#detail/COPERNICUS%2FS2
|
||||
# 10-30m, 10 day repeat
|
||||
bands_s2 = [
|
||||
'B1', # Aerosols
|
||||
'B2', # B
|
||||
'B3', # G
|
||||
'B4', # R
|
||||
'B5', # Red Edge 1 705nm
|
||||
'B6', # Red Edge 2
|
||||
'B7', # Red Edge 3
|
||||
'B8', # NIR
|
||||
'B8A',
|
||||
'B9',
|
||||
'B10',
|
||||
'B11',
|
||||
'B12', # 2190 nm
|
||||
# 'QA10', # empty
|
||||
# 'QA20', # empty
|
||||
'QA60', # cloud
|
||||
]
|
||||
|
||||
# Repeats 6 days?
|
||||
bands_s1 = [
|
||||
'VV',
|
||||
'HH',
|
||||
'VH',
|
||||
'HV',
|
||||
'angle'
|
||||
]
|
||||
|
||||
# https://explorer.earthengine.google.com/#detail/LANDSAT%2FLE7_L1T
|
||||
# 16-60m resolution 16 day repeat
|
||||
bands_l7 = [
|
||||
'B1',
|
||||
'B2',
|
||||
'B3',
|
||||
'B4',
|
||||
'B5',
|
||||
'B6_VCID_1',
|
||||
'B6_VCID_2',
|
||||
'B7',
|
||||
'B8'
|
||||
]
|
||||
|
||||
# https://explorer.earthengine.google.com/#detail/LANDSAT%2FLC8_L1T
|
||||
# 15-100m, 16 day repeat
|
||||
bands_l8 = [
|
||||
'B1',
|
||||
'B2',
|
||||
'B3',
|
||||
'B4',
|
||||
'B5',
|
||||
'B6',
|
||||
'B7',
|
||||
'B8',
|
||||
'B9',
|
||||
'B10',
|
||||
'B11',
|
||||
'BQA'
|
||||
]
|
||||
|
||||
|
||||
def display_ee(geom):
|
||||
"""show earth-engine object in folium map"""
|
||||
import folium
|
||||
# example https://github.com/python-visualization/folium/blob/master/examples/Geopandas.ipynb
|
||||
geojson = geom.getInfo()
|
||||
center = geom.centroid(maxError=1).coordinates().getInfo()
|
||||
|
||||
# folium.initialize_notebook()
|
||||
|
||||
map_osm = folium.Map(location=[center[1], center[0]], tiles='Stamen Terrain')
|
||||
folium.GeoJson(geojson).add_to(map_osm)
|
||||
return map_osm
|
||||
|
||||
# def eethumb(image):
|
||||
# """Show ee image thumbnail in jupyter notebook"""
|
||||
# from IPython.display import HTML
|
||||
# return HTML('<img src="'+image.getThumbUrl()+'"/>')
|
||||
|
||||
|
||||
# A tqdm progress bar for urlretrieve see https://github.com/tqdm/tqdm#hooks-and-callbacks
|
||||
def my_hook(t):
|
||||
"""
|
||||
Wraps tqdm instance. Don't forget to close() or __exit__()
|
||||
the tqdm instance once you're done with it (easiest using `with` syntax).
|
||||
|
||||
Example
|
||||
-------
|
||||
|
||||
>>> with tqdm(...) as t:
|
||||
... reporthook = my_hook(t)
|
||||
... urllib.urlretrieve(..., reporthook=reporthook)
|
||||
|
||||
"""
|
||||
last_b = [0]
|
||||
|
||||
def inner(b=1, bsize=1, tsize=None):
|
||||
"""
|
||||
b : int, optional
|
||||
Number of blocks just transferred [default: 1].
|
||||
bsize : int, optional
|
||||
Size of each block (in tqdm units) [default: 1].
|
||||
tsize : int, optional
|
||||
Total size (in tqdm units). If [default: None] remains unchanged.
|
||||
"""
|
||||
if tsize is not None:
|
||||
t.total = tsize
|
||||
t.update((b - last_b[0]) * bsize)
|
||||
last_b[0] = b
|
||||
|
||||
return inner
|
||||
|
||||
|
||||
# https://github.com/google/earthengine-api/blob/master/python/examples/Image/download.py
|
||||
def download_image(clipped_image, scale=10, crs=4326, name=None, cache_dir=temp_dir, progress_bar=False, report=False):
|
||||
"""Download image from google earth engine"""
|
||||
|
||||
# TODO add progress bar like https://github.com/fchollet/keras/blob/master/keras/utils/data_utils.py#L103
|
||||
path = clipped_image.getDownloadURL({
|
||||
'scale': scale,
|
||||
'crs': 'EPSG:%s' % crs
|
||||
})
|
||||
if name is None:
|
||||
name = clipped_image.getMapId()['mapid']
|
||||
|
||||
# TODO check that it's not too big, but we don't get given clipped size
|
||||
|
||||
filename = '{name:}_{crs:}_{scale:}'.format(
|
||||
name=name, scale=scale, crs=crs)
|
||||
zip_dwn_file = temp_dir.joinpath(filename + '.zip')
|
||||
if report:
|
||||
with tqdm(
|
||||
unit='B', unit_scale=True, miniters=1, mininterval=1,
|
||||
desc=path.split('/')[-1][:40]) as t: # all optional kwargs
|
||||
zip_dwn_file, r = urllib.request.urlretrieve(
|
||||
path, zip_dwn_file, reporthook=my_hook(t))
|
||||
else:
|
||||
zip_dwn_file, r = urllib.request.urlretrieve(path, zip_dwn_file)
|
||||
zip_dwn_file
|
||||
|
||||
# extract
|
||||
zfile = zipfile.ZipFile(zip_dwn_file)
|
||||
extract_dir = cache_dir.joinpath(filename)
|
||||
zfile.extractall(extract_dir)
|
||||
# logger.debug(extract_dir, zip_dwn_file)
|
||||
|
||||
files = [str(f.relpath(extract_dir)) for f in extract_dir.listdir()]
|
||||
logger.debug('Extracted files %s to %s', files, str(extract_dir))
|
||||
return extract_dir, files
|
||||
|
||||
import pyproj
|
||||
import numpy as np
|
||||
|
||||
|
||||
def get_boundary(leak, distance=100, maxError=None):
|
||||
"""get rectangular around geopandas point"""
|
||||
# coords = np.array(leak.geometry.values[0].xy)[:, 0].tolist()
|
||||
# geom = ee.Geometry.Point(coords)
|
||||
# boundary = ee.Geometry.buffer(geometry=geom, distance=distance, maxError=maxError)
|
||||
# rect = boundary.bounds()
|
||||
|
||||
# Here we make a boundary in wgs84 that will make an exact rectangle in
|
||||
# epsg 3857 (aux sphere) so we will get a rectangle when we clip the image
|
||||
point_aux = np.array(leak.geometry.to_crs(epsg=3857).values[0].xy)[:, 0]
|
||||
xMin, yMin, xMax, yMax = point_aux[0] - distance, point_aux[1] - distance, point_aux[0] + distance, point_aux[1] + distance
|
||||
|
||||
# convert to wgs8
|
||||
p0 = pyproj.Proj(init='epsg:%s' % 4326)
|
||||
p1 = pyproj.Proj(init='epsg:%s' % 3857)
|
||||
bound_grid2 = pyproj.transform(p1, p0, [xMin, xMax], [yMin, yMax])
|
||||
bound_grid2 = np.array(bound_grid2).T
|
||||
[xMin, yMin], [xMax, yMax] = bound_grid2.min(0), bound_grid2.max(0)
|
||||
|
||||
# make into earth engine rec
|
||||
rect = ee.Geometry.Rectangle([xMin, yMin, xMax, yMax])
|
||||
return rect
|
||||
|
||||
# boundary = get_boundary(leak.geometry)
|
||||
# # TEST TODO
|
||||
# b = boundary.getInfo()
|
||||
# import pyproj
|
||||
# p0=pyproj.Proj(init='epsg:%s'%4326)
|
||||
# p1=pyproj.Proj(init='epsg:%s'%3857)
|
||||
# bb=np.array(b['coordinates'][0])
|
||||
# bound_grid = pyproj.transform(p0,p1,bb[:,0],bb[:,1])
|
||||
# # get distance in aux grid
|
||||
# bg2=np.array(bound_grid).T
|
||||
# assert bg2.max(0)-bg2.min(0)==[distance, distance]
|
||||
# assert (bg2.max(0)-bg2.min(0))/min_resolution==pixel_length
|
||||
#
|
||||
# def image2array(image, point, crs=4326, resolution_min=resolution_min, pixel_length=pixel_length, name=None):
|
||||
#
|
||||
# scale = resolution_min
|
||||
#
|
||||
# # get image
|
||||
# path,files=download_image(image, scale=scale, crs=crs, name=name)
|
||||
#
|
||||
# data = tifs2np(path,files,bands=bands)
|
||||
#
|
||||
# # now if the size is wrong let's interp it
|
||||
# if data.shape[-2]!=pixel_length or data.shape[-1]!=pixel_length:
|
||||
# data = np.array([sp.misc.imresize(x,size=(pixel_length,pixel_length),interp='cubic', mode='F') for x in data])
|
||||
# return data
|
||||
# test
|
||||
# data = tifs2np(path,files)
|
||||
|
||||
import scipy as sp
|
||||
|
||||
|
||||
def tifs2np(path, files, pixel_length=None, bands=bands_s2):
|
||||
"""Convert tifs to numpy array"""
|
||||
tifs = [f for f in files if f.endswith('.tif')]
|
||||
|
||||
if pixel_length:
|
||||
pixel_length = int(pixel_length)
|
||||
|
||||
channels = {}
|
||||
for tif in tifs:
|
||||
band = tif.split('.')[-2]
|
||||
# read tif as float32
|
||||
x = imread(path.joinpath(tif), mode='F')
|
||||
if pixel_length:
|
||||
if x.shape[-2] != pixel_length or x.shape[-1] != pixel_length:
|
||||
logger.warn('warning had to reshape band %s from %s to %s' % (band, x.shape, pixel_length))
|
||||
x = sp.misc.imresize(x, size=(pixel_length, pixel_length), interp='cubic', mode='F')
|
||||
channels[band] = x
|
||||
|
||||
if not pixel_length:
|
||||
pixel_length = x.shape[1]
|
||||
|
||||
logger.debug('keys %s', (channels.keys()))
|
||||
data = []
|
||||
for band in bands:
|
||||
if band not in channels:
|
||||
channels[band] = np.zeros((pixel_length, pixel_length))
|
||||
data.append(channels[band])
|
||||
return np.array(data)
|
||||
@@ -0,0 +1,85 @@
|
||||
import pylab
|
||||
import shapely
|
||||
import numpy as np
|
||||
|
||||
DEFAULT_MIN_SIZE = 10
|
||||
DEFAULT_SAMPLE = 10
|
||||
|
||||
|
||||
def diffxy(line):
|
||||
"""diff of x,y array in terms of dist."""
|
||||
# check it the right way round, otherwise I'll get false results
|
||||
if len(line.shape) == 1:
|
||||
# return np.sqrt(np.square(line[0])+np.square(line[1]))
|
||||
line = line.reshape(1, -1)
|
||||
|
||||
r, c = line.shape
|
||||
if r < c:
|
||||
line = line.T
|
||||
|
||||
dists = pylab.distances_along_curve(line)
|
||||
dists = np.hstack((0, dists))
|
||||
return dists
|
||||
|
||||
# note this is about 100x as fast as shapely
|
||||
|
||||
|
||||
def resample_polygon(line, n=np.floor(DEFAULT_MIN_SIZE / DEFAULT_SAMPLE), records=[], step=True):
|
||||
"""
|
||||
Interp based on cumdist.
|
||||
|
||||
we loose one point of the length but add it again
|
||||
|
||||
n is samples distance,unless step=False when it's number of points
|
||||
|
||||
good documentation about the errors in this method:
|
||||
http://stackoverflow.com/questions/4052225/how-to-equidistant-resample-a-line-or-curve.
|
||||
"""
|
||||
l = np.cumsum(diffxy(line))
|
||||
|
||||
if step:
|
||||
if l[-1] < n:
|
||||
return line, records
|
||||
else:
|
||||
if len(line) < n:
|
||||
return line, records
|
||||
|
||||
x = line[:, 0]
|
||||
y = line[:, 1]
|
||||
if step:
|
||||
ll = np.arange(min(l), max(l), n)
|
||||
else:
|
||||
ll = np.linspace(min(l), max(l), n)
|
||||
|
||||
# linear piecewise
|
||||
xx = np.interp(ll, l, x)
|
||||
yy = np.interp(ll, l, y)
|
||||
|
||||
# from scipy.interpolate import interp1d
|
||||
# kind='cubic' # cubic
|
||||
# fx=interp1d(l,x,kind=kind)
|
||||
# xx=fx(ll)
|
||||
# fy=interp1d(l,y,kind=kind)
|
||||
# yy=fy(ll)
|
||||
# diffxy(ar(zip(xx,yy)))[1:].std()
|
||||
# linear: 2016, nearest=2028, slinear=2016, quad: 1E20, cubic
|
||||
|
||||
# from scipy.interpolate import UnivariateSpline
|
||||
|
||||
# now add original end back on so as not too lose data
|
||||
xx[-1] = x[-1]
|
||||
yy[-1] = y[-1]
|
||||
|
||||
# plot(xx,yy);plot(x,y)
|
||||
line2 = np.array(list(zip(xx, yy)))
|
||||
if records != []:
|
||||
rrecords = []
|
||||
for i in xrange(records.shape[1]):
|
||||
r = records[:, i]
|
||||
rr = np.interp(ll, l, r)
|
||||
# rr=np.append(r[0],rr)
|
||||
rrecords.append(rr)
|
||||
rrecords = np.dstack(rrecords)[0]
|
||||
return line2, rrecords
|
||||
else:
|
||||
return line2
|
||||
@@ -0,0 +1,3 @@
|
||||
from .dice_loss import *
|
||||
from .filters import *
|
||||
from .multichannel_imagedatagenerator import *
|
||||
@@ -0,0 +1,20 @@
|
||||
"""
|
||||
Here is a dice loss for keras which is smoothed to approximate a linear (L1) loss.
|
||||
It ranges from 1 to 0 (no error), and returns results similar to binary crossentropy
|
||||
"""
|
||||
|
||||
# define custom loss and metric functions
|
||||
|
||||
from keras import backend as K
|
||||
|
||||
def dice_coef(y_true, y_pred, smooth=1):
|
||||
"""
|
||||
Dice = (2*|X & Y|)/ (|X|+ |Y|)
|
||||
= 2*sum(|A*B|)/(sum(A^2)+sum(B^2))
|
||||
ref: https://arxiv.org/pdf/1606.04797v1.pdf
|
||||
"""
|
||||
intersection = K.sum(K.abs(y_true * y_pred), axis=-1)
|
||||
return (2. * intersection + smooth) / (K.sum(K.square(y_true),-1) + K.sum(K.square(y_pred),-1) + smooth)
|
||||
|
||||
def dice_coef_loss(y_true, y_pred):
|
||||
return 1-dice_coef(y_true, y_pred)
|
||||
@@ -0,0 +1,167 @@
|
||||
"""
|
||||
Filters for input data
|
||||
"""
|
||||
import numpy as np
|
||||
import arrow
|
||||
import logging
|
||||
logger = logging.getLogger('leak_helpers.modelling.filters')
|
||||
|
||||
|
||||
def normalise_bands(X):
|
||||
"""Each band becomes -1 to 1 and 0 centered"""
|
||||
for i in range(X.shape[1]):
|
||||
|
||||
x = X[:, i, :, :]
|
||||
if x.std() > 0:
|
||||
X[:, i, :, :] = (x - x.mean()) / (x.max() - x.min())
|
||||
return X
|
||||
|
||||
|
||||
def is_not_cloudy(md, max_cover=0.3):
|
||||
if 'CLOUD_COVER' in md['image']['properties']:
|
||||
return md['image']['properties']['CLOUD_COVER'] / 100.0 < max_cover
|
||||
elif 'CLOUDY_PIXEL_PERCENTAGE' in md['image']['properties']:
|
||||
return md['image']['properties']['CLOUDY_PIXEL_PERCENTAGE'] / 100.0 < max_cover
|
||||
else:
|
||||
return True
|
||||
|
||||
|
||||
def is_not_center_cloudy(X):
|
||||
"""
|
||||
Check if the center pixel have the cloud mask on them.
|
||||
|
||||
It checks the last band which is cloud in sentinel-2 and cloud/snow/etc in landsat-8.
|
||||
"""
|
||||
# TODO for landsat-8 I need to look at std
|
||||
cut = int(np.ceil(X.shape[-2] / 2.0 - 1)) # if it's 25 wide cut at 12, if it's 24, at 11
|
||||
return X[:, -1, cut:-cut, cut:-cut].reshape((X.shape[0], -1)).any(-1) == False
|
||||
|
||||
|
||||
def is_image_within(md, seconds=60 * 60 * 24):
|
||||
"""Check image is within X time."""
|
||||
t_image = arrow.get(md['image']['properties']['system:time_end'] / 1000)
|
||||
t_leak = arrow.get(md['leak']['features'][0]['properties']['REPO_Date'])
|
||||
seconds_before_leak = (t_leak - t_image).total_seconds()
|
||||
return seconds_before_leak < seconds
|
||||
|
||||
|
||||
def hash_rows(X_train):
|
||||
return [hash(X_train[i].tobytes()) for i in range(len(X_train))]
|
||||
|
||||
|
||||
def is_not_dup(X1):
|
||||
"""False for the duplicates (True for the first instance)"""
|
||||
n = hash_rows(X1)
|
||||
return np.array([(n[i] not in n[:i]) for i in range(len(n))])
|
||||
|
||||
|
||||
def is_leak(md):
|
||||
"""Check is says LEAK in one of the text fields"""
|
||||
props = md['leak']['features'][0]['properties']
|
||||
if 'WA' not in props['leak_id']:
|
||||
return True
|
||||
else:
|
||||
return np.any(["LEAK" in str(v) for v in md['leak']['features'][0]['properties'].values()])
|
||||
|
||||
|
||||
from sklearn.model_selection import train_test_split
|
||||
import pandas as pd
|
||||
from sklearn.utils import shuffle
|
||||
|
||||
|
||||
def filter_data(X_raw, y_raw, metadatas, max_cloud_cover=1, timespan_before=np.inf, random_seed=0, normalized=True, balanced_classes=True, filter_center_cloudy=False):
|
||||
# filter based on cloud and timespan
|
||||
filtr = np.array([
|
||||
[
|
||||
is_not_cloudy(metadata, max_cover=max_cloud_cover),
|
||||
is_leak(metadata),
|
||||
is_image_within(metadata, timespan_before),
|
||||
] for metadata in metadatas
|
||||
])
|
||||
if filter_center_cloudy:
|
||||
incc = is_not_center_cloudy(X_raw)
|
||||
else:
|
||||
incc = np.ones((len(X_raw)))
|
||||
filtr = np.hstack([
|
||||
filtr,
|
||||
incc.reshape((-1, 1)),
|
||||
is_not_dup(X_raw).reshape((-1, 1))
|
||||
])
|
||||
filt = filtr.all(-1)
|
||||
|
||||
# # QC filter
|
||||
# df_filt = pd.DataFrame(filtr, columns=['is_not_cloudy','is_image_within','is_not_center_cloudy'])
|
||||
# print('filter', df_filt.sum(0), len(df_filt))
|
||||
# print('"before" images passing time filter', df_filt.is_image_within.sum()-len(df_filt)/2)
|
||||
|
||||
metadata_filtered = [metadatas[i] for i in range(len(filt)) if filt[i]]
|
||||
X = X_raw[filt]
|
||||
y = y_raw[filt]
|
||||
|
||||
# now because I want balanced data I have to filter them by id now
|
||||
if balanced_classes:
|
||||
|
||||
allowed_ids = set([metadata_filtered[i]['leak']['features'][0]['properties']['leak_id'] for i in range(len(y)) if y[i]])
|
||||
logger.debug('allowed_ids', len(allowed_ids))
|
||||
|
||||
# there are move of class 0, so we mask them until we reach balance
|
||||
mask = np.ones((len(y)), dtype=np.bool)
|
||||
for i in range(len(mask)):
|
||||
if y[mask].mean() >= 0.5:
|
||||
# stop because it's balanced
|
||||
break
|
||||
if not y[i]:
|
||||
leak_id = metadata_filtered[i]['leak']['features'][0]['properties']['leak_id']
|
||||
# we target the onces that pair with the removed class 1's
|
||||
if leak_id not in allowed_ids:
|
||||
mask[i] = False
|
||||
|
||||
metadata_filtered = [metadata_filtered[i] for i in range(len(mask)) if mask[i]]
|
||||
X = X[mask]
|
||||
y = y[mask]
|
||||
X.shape, y.shape, len(metadata_filtered)
|
||||
|
||||
if np.abs(0.5 - y.mean()) > 0.1:
|
||||
logger.error('balanced classes should have an even number of each')
|
||||
|
||||
if normalized:
|
||||
X = normalise_bands(X)
|
||||
|
||||
# shuffle
|
||||
return X, y, metadata_filtered
|
||||
|
||||
|
||||
def filter_split_data(X_raw, y_raw, metadatas, max_cloud_cover=1, timespan_before=np.inf, test_fraction=0.3, val_fraction=0.3, random_seed=0, normalized=True, balanced_classes=True, filter_center_cloudy=False):
|
||||
X, y, metadata_filtered = filter_data(X_raw, y_raw, metadatas, max_cloud_cover=max_cloud_cover, timespan_before=timespan_before, random_seed=random_seed, normalized=normalized, balanced_classes=balanced_classes, filter_center_cloudy=filter_center_cloudy)
|
||||
|
||||
X, y, metadata_filtered=shuffle(X, y, metadata_filtered, random_state=random_seed)
|
||||
|
||||
X_train, X_test, y_train, y_test, metadata_train, metadata_test=train_test_split(
|
||||
X, y, metadata_filtered, test_size=test_fraction, random_state=random_seed)
|
||||
|
||||
X_train, X_val, y_train, y_val, metadata_train, metadata_val=train_test_split(
|
||||
X_train, y_train, metadata_train, test_size=val_fraction, random_state=random_seed)
|
||||
# print(X_train.shape,y_train.shape, len(metadata_train))
|
||||
# print(X_test.shape,y_test.shape, len(metadata_test))
|
||||
# print(X_val.shape,y_val.shape, len(metadata_val))
|
||||
|
||||
return X_train, y_train, metadata_train, X_val, y_val, metadata_val, X_test, y_test, metadata_test
|
||||
|
||||
# X_raw = np.random.random((100,10,10,10))
|
||||
# y_raw = np.random.random(100)>0.5
|
||||
# metadata = [dict(leak_id=i,REPO_Date=arrow.get().timestamp) for i in range(100)]
|
||||
# X_train, y_train, metadata_train, X_val, y_val, metadata_val, X_test, y_test, metadata_test = filter_split_data(
|
||||
# X_raw,
|
||||
# y_raw,
|
||||
# metadatas,
|
||||
# max_cloud_cover=0.3,
|
||||
# timespan_before=60*60*24*3,
|
||||
# test_fraction=0.3,
|
||||
# random_seed=0,
|
||||
# balanced_classes=True,
|
||||
# normalized=False,
|
||||
# filter_center_cloudy=True,
|
||||
# )
|
||||
# print(X_train.shape,y_train.shape, len(metadata_train))
|
||||
# print(X_test.shape,y_test.shape, len(metadata_test))
|
||||
# print(X_val.shape,y_val.shape, len(metadata_val))
|
||||
@@ -0,0 +1,152 @@
|
||||
|
||||
# mod to Allow the ImageDateGenerator to have multiple channels instead of just 1,3,4
|
||||
|
||||
# modified from https://github.com/fchollet/keras/blob/master/keras/preprocessing/image.py
|
||||
from keras.preprocessing.image import Iterator
|
||||
from keras import backend as K
|
||||
from keras.preprocessing.image import ImageDataGenerator as _ImageDataGenerator
|
||||
from path import Path
|
||||
from scipy import linalg
|
||||
import numpy as np
|
||||
|
||||
class NumpyArrayIterator(Iterator):
|
||||
|
||||
def __init__(self, x, y, image_data_generator,
|
||||
batch_size=32, shuffle=False, seed=None,
|
||||
dim_ordering='default',
|
||||
save_to_dir=None, save_prefix='', save_format='jpeg'):
|
||||
if y is not None and len(x) != len(y):
|
||||
raise ValueError('X (images tensor) and y (labels) '
|
||||
'should have the same length. '
|
||||
'Found: X.shape = %s, y.shape = %s' %
|
||||
(np.asarray(x).shape, np.asarray(y).shape))
|
||||
if dim_ordering == 'default':
|
||||
dim_ordering = K.image_dim_ordering()
|
||||
self.x = np.asarray(x)
|
||||
if self.x.ndim != 4:
|
||||
raise ValueError('Input data in `NumpyArrayIterator` '
|
||||
'should have rank 4. You passed an array '
|
||||
'with shape', self.x.shape)
|
||||
channels_axis = 3 if dim_ordering == 'tf' else 1
|
||||
# if self.x.shape[channels_axis] not in {1, 3, 4}:
|
||||
# raise ValueError('NumpyArrayIterator is set to use the '
|
||||
# 'dimension ordering convention "' + dim_ordering + '" '
|
||||
# '(channels on axis ' + str(channels_axis) + '), i.e. expected '
|
||||
# 'either 1, 3 or 4 channels on axis ' + str(channels_axis) + '. '
|
||||
# 'However, it was passed an array with shape ' + str(self.x.shape) +
|
||||
# ' (' + str(self.x.shape[channels_axis]) + ' channels).')
|
||||
if y is not None:
|
||||
self.y = np.asarray(y)
|
||||
else:
|
||||
self.y = None
|
||||
self.image_data_generator = image_data_generator
|
||||
self.dim_ordering = dim_ordering
|
||||
self.save_to_dir = save_to_dir
|
||||
self.save_prefix = save_prefix
|
||||
self.save_format = save_format
|
||||
super(NumpyArrayIterator, self).__init__(x.shape[0], batch_size, shuffle, seed)
|
||||
|
||||
def next(self):
|
||||
# for python 2.x.
|
||||
# Keeps under lock only the mechanism which advances
|
||||
# the indexing of each batch
|
||||
# see http://anandology.com/blog/using-iterators-and-generators/
|
||||
with self.lock:
|
||||
index_array, current_index, current_batch_size = next(self.index_generator)
|
||||
# The transformation of images is not under thread lock
|
||||
# so it can be done in parallel
|
||||
batch_x = np.zeros(tuple([current_batch_size] + list(self.x.shape)[1:]))
|
||||
for i, j in enumerate(index_array):
|
||||
x = self.x[j]
|
||||
x = self.image_data_generator.random_transform(x.astype('float32'))
|
||||
x = self.image_data_generator.standardize(x)
|
||||
batch_x[i] = x
|
||||
if self.save_to_dir:
|
||||
for i in range(current_batch_size):
|
||||
img = array_to_img(batch_x[i], self.dim_ordering, scale=True)
|
||||
fname = '{prefix}_{index}_{hash}.{format}'.format(prefix=self.save_prefix,
|
||||
index=current_index + i,
|
||||
hash=np.random.randint(1e4),
|
||||
format=self.save_format)
|
||||
img.save(os.path.join(self.save_to_dir, fname))
|
||||
if self.y is None:
|
||||
return batch_x
|
||||
batch_y = self.y[index_array]
|
||||
return batch_x, batch_y
|
||||
|
||||
|
||||
class ImageDataGenerator(_ImageDataGenerator):
|
||||
|
||||
def flow(self, X, y=None, batch_size=32, shuffle=True, seed=None,
|
||||
save_to_dir=None, save_prefix='', save_format='jpeg'):
|
||||
return NumpyArrayIterator(
|
||||
X, y, self,
|
||||
batch_size=batch_size,
|
||||
shuffle=shuffle,
|
||||
seed=seed,
|
||||
dim_ordering=self.dim_ordering,
|
||||
save_to_dir=save_to_dir,
|
||||
save_prefix=save_prefix,
|
||||
save_format=save_format)
|
||||
|
||||
def fit(self, x,
|
||||
augment=False,
|
||||
rounds=1,
|
||||
seed=None):
|
||||
"""Required for featurewise_center, featurewise_std_normalization
|
||||
and zca_whitening.
|
||||
# Arguments
|
||||
x: Numpy array, the data to fit on. Should have rank 4.
|
||||
In case of grayscale data,
|
||||
the channels axis should have value 1, and in case
|
||||
of RGB data, it should have value 3.
|
||||
augment: Whether to fit on randomly augmented samples
|
||||
rounds: If `augment`,
|
||||
how many augmentation passes to do over the data
|
||||
seed: random seed.
|
||||
# Raises
|
||||
ValueError: in case of invalid input `x`.
|
||||
"""
|
||||
x = np.asarray(x)
|
||||
if x.ndim != 4:
|
||||
raise ValueError('Input to `.fit()` should have rank 4. '
|
||||
'Got array with shape: ' + str(x.shape))
|
||||
# if x.shape[self.channel_axis] not in {1, 3, 4}:
|
||||
# raise ValueError(
|
||||
# 'Expected input to be images (as Numpy array) '
|
||||
# 'following the dimension ordering convention "' + self.dim_ordering + '" '
|
||||
# '(channels on axis ' + str(self.channel_axis) + '), i.e. expected '
|
||||
# 'either 1, 3 or 4 channels on axis ' + str(self.channel_axis) + '. '
|
||||
# 'However, it was passed an array with shape ' + str(x.shape) +
|
||||
# ' (' + str(x.shape[self.channel_axis]) + ' channels).')
|
||||
|
||||
if seed is not None:
|
||||
np.random.seed(seed)
|
||||
|
||||
x = np.copy(x)
|
||||
if augment:
|
||||
ax = np.zeros(tuple([rounds * x.shape[0]] + list(x.shape)[1:]))
|
||||
for r in range(rounds):
|
||||
for i in range(x.shape[0]):
|
||||
ax[i + r * x.shape[0]] = self.random_transform(x[i])
|
||||
x = ax
|
||||
|
||||
if self.featurewise_center:
|
||||
self.mean = np.mean(x, axis=(0, self.row_axis, self.col_axis))
|
||||
broadcast_shape = [1, 1, 1]
|
||||
broadcast_shape[self.channel_axis - 1] = x.shape[self.channel_axis]
|
||||
self.mean = np.reshape(self.mean, broadcast_shape)
|
||||
x -= self.mean
|
||||
|
||||
if self.featurewise_std_normalization:
|
||||
self.std = np.std(x, axis=(0, self.row_axis, self.col_axis))
|
||||
broadcast_shape = [1, 1, 1]
|
||||
broadcast_shape[self.channel_axis - 1] = x.shape[self.channel_axis]
|
||||
self.std = np.reshape(self.std, broadcast_shape)
|
||||
x /= (self.std + K.epsilon())
|
||||
|
||||
if self.zca_whitening:
|
||||
flat_x = np.reshape(x, (x.shape[0], x.shape[1] * x.shape[2] * x.shape[3]))
|
||||
sigma = np.dot(flat_x.T, flat_x) / flat_x.shape[0]
|
||||
u, s, _ = linalg.svd(sigma)
|
||||
self.principal_components = np.dot(np.dot(u, np.diag(1. / np.sqrt(s + 10e-7))), u.T)
|
||||
@@ -0,0 +1,31 @@
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
|
||||
# TODO show band labels?
|
||||
def imshow_bands(X_img, figsize=(14, 14), interpolation=None, labels=None):
|
||||
"""show mosiac of input bands."""
|
||||
plt.figure(figsize=figsize)
|
||||
|
||||
f, axarr = plt.subplots(4, 4)
|
||||
for j in range(X_img.shape[0]):
|
||||
plt.subplot(4, 4, j + 1)
|
||||
# plt.subplot(j)
|
||||
# ax=axarr.flatten()[i]
|
||||
plt.imshow(X_img[j, :, :], interpolation=interpolation)
|
||||
plt.grid(False)
|
||||
plt.gca().set_axis_off()
|
||||
if labels:
|
||||
label = labels[j]
|
||||
else:
|
||||
label = ''
|
||||
plt.title('#%i - %s' % (j + 1, label))
|
||||
# plt.show()
|
||||
|
||||
plt.subplot(4, 4, j + 2)
|
||||
plt.grid(False)
|
||||
plt.gca().set_axis_off()
|
||||
plt.subplot(4, 4, j + 3)
|
||||
plt.grid(False)
|
||||
plt.gca().set_axis_off()
|
||||
plt.suptitle('Image bands')
|
||||
return plt.gca()
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,920 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This notebook scrapes satellite images for each leak repair. For each location it gets a NxM rectangle around the leak before and after it was repaired. Then it collated all the data into h5 files and all the metadata into json files.\n",
|
||||
"\n",
|
||||
"It takes days to run because of rate limiting on the google earth api. Because of limited satelite coverage you might find matches for only 10% of the leaks.\n",
|
||||
"\n",
|
||||
"## Modifying\n",
|
||||
"\n",
|
||||
"- make sure google earth is setup\n",
|
||||
"- load leaks, so they pass the asserts\n",
|
||||
"- change params\n",
|
||||
"- run rest of cells"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:09.472506Z",
|
||||
"start_time": "2017-03-18T18:08:07.757873+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"'%.4f'"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from path import Path\n",
|
||||
"import arrow\n",
|
||||
"import json\n",
|
||||
"import pytz\n",
|
||||
"from pprint import pprint\n",
|
||||
"from tqdm import tqdm_notebook as tqdm\n",
|
||||
"import re, os, collections, itertools, uuid, logging\n",
|
||||
"import tempfile\n",
|
||||
"import tables\n",
|
||||
"\n",
|
||||
"import zipfile\n",
|
||||
"import urllib\n",
|
||||
"\n",
|
||||
"import ee\n",
|
||||
"import pyproj\n",
|
||||
"import numpy as np\n",
|
||||
"import scipy as sp\n",
|
||||
"import pandas as pd\n",
|
||||
"import geopandas as gpd\n",
|
||||
"from matplotlib import pyplot as plt\n",
|
||||
"import seaborn as sns\n",
|
||||
"\n",
|
||||
"plt.rcParams['figure.figsize'] = (15, 5) # bigger plots\n",
|
||||
"plt.style.use('fivethirtyeight')\n",
|
||||
"%matplotlib inline\n",
|
||||
"%precision 4"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:09.477069Z",
|
||||
"start_time": "2017-03-18T18:08:09.474514+08:00"
|
||||
},
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# %load_ext autoreload\n",
|
||||
"# %autoreload 2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:12.811993Z",
|
||||
"start_time": "2017-03-18T18:08:09.479209+08:00"
|
||||
},
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"helper_dir = str(Path('..').abspath())\n",
|
||||
"if helper_dir not in os.sys.path:\n",
|
||||
" os.sys.path.append(helper_dir)\n",
|
||||
" \n",
|
||||
"from leak_helpers.earth_engine import display_ee, get_boundary, tifs2np, bands_s2, download_image, bands_s2, bands_s1, bands_l7, bands_l8"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:12.848915Z",
|
||||
"start_time": "2017-03-18T18:08:12.813510+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(Path('/tmp/testing_earth_engine-l7-AUTX-48wtrw24-20170318-10-08-12'),\n",
|
||||
" Path('../data/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2'),\n",
|
||||
" Path('../data/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2/ee_l7_AUTX-leaks_cache_v2'))"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"\n",
|
||||
"crs_grid = 3857\n",
|
||||
"notebook_name='testing_earth_engine-l7-AUTX'\n",
|
||||
"ts=arrow.utcnow().format('YYYYMMDD-HH-mm-ss')\n",
|
||||
"data_dir = Path('../data/')\n",
|
||||
"bands = bands_l7\n",
|
||||
"\n",
|
||||
"# since the lowest res band is 60m and I want to capture neighbours I should get 6+ pixels\n",
|
||||
"pixel_length = 25.0\n",
|
||||
"resolution_min = 15.0 # m\n",
|
||||
"time_bin_delta = 60*60*24*28 # how long before a leak to look (in seconds)\n",
|
||||
"# TODO get closest but let me filter for time\n",
|
||||
"\n",
|
||||
"# init\n",
|
||||
"temp_dir = Path(tempfile.mkdtemp(prefix=notebook_name+'-', suffix='-'+ts))\n",
|
||||
"# output_dir = data_dir.joinpath('{ts:}_{notebook_name:}'.format(ts=ts,notebook_name=notebook_name))\n",
|
||||
"output_dir = Path('../../data/scraped_satellite_images/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2')\n",
|
||||
"cache_dir = output_dir.joinpath('ee_l7_AUTX-leaks_cache_v2')\n",
|
||||
"\n",
|
||||
"output_dir.makedirs_p()\n",
|
||||
"temp_dir.makedirs_p()\n",
|
||||
"cache_dir.makedirs_p()\n",
|
||||
"\n",
|
||||
"logger = logging.getLogger(notebook_name)\n",
|
||||
"logger.setLevel(logging.WARN)\n",
|
||||
"\n",
|
||||
"crs_grid_proj = pyproj.Proj('+init=epsg:%s'%crs_grid)\n",
|
||||
"\n",
|
||||
"temp_dir, output_dir, cache_dir"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:12.870492Z",
|
||||
"start_time": "2017-03-18T18:08:12.850253+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"metadata_file = output_dir.joinpath('script_metadata.json')\n",
|
||||
"\n",
|
||||
"# write metadata to json\n",
|
||||
"metadata = dict(\n",
|
||||
" pixel_length=pixel_length,\n",
|
||||
" resolution_min=resolution_min,\n",
|
||||
" bands=bands,\n",
|
||||
" ts=ts,\n",
|
||||
" notebook_name=notebook_name,\n",
|
||||
" crs_grid=crs_grid,\n",
|
||||
" cache_dir=str(cache_dir),\n",
|
||||
" temp_dir=str(temp_dir),\n",
|
||||
" output_dir=str(output_dir),\n",
|
||||
")\n",
|
||||
"json.dump(metadata, open(metadata_file,'w'))\n",
|
||||
"\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# earth engine\n",
|
||||
"\n",
|
||||
"Setup instructions here\n",
|
||||
"- first need to apply for an account and wait ~ 1day\n",
|
||||
"- https://developers.google.com/earth-engine/python_install#setting-up-authentication-credentials\n",
|
||||
"\n",
|
||||
"Refs:\n",
|
||||
"- api https://developers.google.com/earth-engine/\n",
|
||||
"- code examples https://code.earthengine.google.com/\n",
|
||||
"- sentinel1 https://developers.google.com/earth-engine/sentinel1\n",
|
||||
" - `ee.ImageCollection('COPERNICUS/S2_GRD');`\n",
|
||||
" - `ee.ImageCollection('COPERNICUS/S1_GRD');`\n",
|
||||
"- keras and google earth https://github.com/patrick-dd/landsat-landstats"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:15.140158Z",
|
||||
"start_time": "2017-03-18T18:08:12.872467+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"ok\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# test earth-engine setup\n",
|
||||
"from oauth2client import crypt # should have not error\n",
|
||||
"import ee\n",
|
||||
"ee.Initialize() # should give no errors, if so follow instructions\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# test\n",
|
||||
"image = ee.Image('srtm90_v4')\n",
|
||||
"assert image.getInfo()=={'type': 'Image', 'properties': {'system:time_start': 950227200000, 'system:asset_size': 18827626666, 'system:time_end': 951177600000}, 'bands': [{'data_type': {'type': 'PixelType', 'max': 32767, 'min': -32768, 'precision': 'int'}, 'crs': 'EPSG:4326', 'id': 'elevation', 'dimensions': [432000, 144000], 'crs_transform': [0.000833333333333, 0.0, -180.0, 0.0, -0.000833333333333, 60.0]}], 'id': 'srtm90_v4', 'version': 1463778555689000}\n",
|
||||
"print('ok')\n",
|
||||
"\n",
|
||||
"# ee.Geometry.Point([117.21079620254062, -30.94712385398404])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Load leaks"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:17.925179Z",
|
||||
"start_time": "2017-03-18T18:08:15.141537+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "NameError",
|
||||
"evalue": "name 'leaks' is not defined",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[0;32m<ipython-input-7-c524ae34d362>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0mleaks_ATX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'leak_id'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mleaks_ATX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mOBJECTID\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m'ATX-%s'\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;31m# leaks=leaks_ATX\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mleaks\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mleaks\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mleak_id\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 12\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mleaks_ATX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;31mNameError\u001b[0m: name 'leaks' is not defined"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# load wa leaks\n",
|
||||
"leaks_ATX = gpd.read_file(data_dir.joinpath('austin_leaks/derived/austin_leaks-repairs.geojson'))\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# they have to be after launch\n",
|
||||
"s3_launch_ts=pd.Timestamp('1 Jan 1999')\n",
|
||||
"leaks_ATX = leaks_ATX[pd.to_datetime(leaks_ATX.COMPDTTM)>=s3_launch_ts]\n",
|
||||
"leaks_ATX['REPO_Date']=leaks_ATX['COMPDTTM']\n",
|
||||
"leaks_ATX['leak_id']=leaks_ATX.OBJECTID.apply(lambda x:'ATX-%s'%x)\n",
|
||||
"# leaks=leaks_ATX\n",
|
||||
"leaks.index = leaks.leak_id\n",
|
||||
"len(leaks_ATX)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:17.925652Z",
|
||||
"start_time": "2017-03-18T10:08:09.292Z"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# choose one leak for now\n",
|
||||
"leak = leaks_ATX.sample()\n",
|
||||
"leak"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-01-15T10:46:31.548298",
|
||||
"start_time": "2017-01-15T10:46:31.546367"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"# Fetching sentinal-1 and sentinel 2 images\n",
|
||||
"\n",
|
||||
"For a leak repair, grab the image before and after it"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Note roughly 10% have results for a 1 day temporal bin"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:28.842346Z",
|
||||
"start_time": "2017-03-18T18:08:28.836767+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def get_cached_ids():\n",
|
||||
" cache_dirs = [str(f.relpath(cache_dir)).split('_')[0] for f in cache_dir.listdir()]\n",
|
||||
" return cache_dirs\n",
|
||||
"\n",
|
||||
"def init_cache(leak_id):\n",
|
||||
" \"\"\"We will cache downloads in folders like 'id_after'\"\"\"\n",
|
||||
" if leak_id:\n",
|
||||
" cache_subdir = cache_dir.joinpath(leak_id+'_after')\n",
|
||||
" cache_subdir.makedirs_p()\n",
|
||||
" cache_subdir = cache_dir.joinpath(leak_id+'_before')\n",
|
||||
" cache_subdir.makedirs_p()\n",
|
||||
" return get_cached_ids()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"For each point\n",
|
||||
"- find the nearest image before the repair\n",
|
||||
"- and the soonest image after repair\n",
|
||||
"- save a part of each with metadata\n",
|
||||
"\n",
|
||||
"Later we can filter, interpolate, and read into numpy arrays"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:30.447941Z",
|
||||
"start_time": "2017-03-18T18:08:30.445152+08:00"
|
||||
},
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"distance = resolution_min*(pixel_length/2.0-0.5)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T10:08:17.926512Z",
|
||||
"start_time": "2017-03-18T10:08:10.530Z"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"\n",
|
||||
"# test with one image\n",
|
||||
"for i in [10,50,100,1000,2000]:\n",
|
||||
" leak=leaks_ATX.iloc[[i]]\n",
|
||||
" leak_id = str(leak.OBJECTID.values[0])\n",
|
||||
"\n",
|
||||
" repo_date_ts = arrow.get(leak.REPO_Date.values[0]).timestamp\n",
|
||||
" boundary = get_boundary(leak, distance=distance)\n",
|
||||
" sentinel2_before = ee.ImageCollection('LANDSAT/LE7_L1T')\\\n",
|
||||
" .filterBounds(boundary)\\\n",
|
||||
" .filterDate(933828614605,1488776737937)\\\n",
|
||||
" .sort('system:time_start', opt_ascending=False) # first will be latest\n",
|
||||
" image = ee.Image(sentinel2_before.first()).clip(boundary)\n",
|
||||
" image.getInfo()\n",
|
||||
" name=leak_id+'_after'\n",
|
||||
" path,files=download_image(\n",
|
||||
" image, \n",
|
||||
" scale=resolution_min, \n",
|
||||
" crs=crs_grid, \n",
|
||||
" name=name,\n",
|
||||
" cache_dir=cache_dir\n",
|
||||
" )\n",
|
||||
" data = tifs2np(path,files,bands=bands_l7)\n",
|
||||
" print(i, [(d.shape,d.sum()) for d in data])\n",
|
||||
" for d in data:\n",
|
||||
" assert d.shape[0]==pixel_length, 'the downloaded image is the wrong size, tweak distance'\n",
|
||||
" assert d.shape[1]==pixel_length\n",
|
||||
" assert np.sum(data)!=0,'should not be empty (make sure you are using the right bands)'"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-15T01:59:52.869742Z",
|
||||
"start_time": "2017-03-15T09:59:52.865652+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-15T14:34:52.122574Z",
|
||||
"start_time": "2017-03-15T22:34:52.120409+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-19T08:08:16.798286Z",
|
||||
"start_time": "2017-03-18T18:08:31.789754+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"application/vnd.jupyter.widget-view+json": {
|
||||
"model_id": "ecbee146c27a476ca62d495430ac5cbe"
|
||||
}
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"17868 Earth Engine memory capacity exceeded.\n",
|
||||
"17869 Earth Engine memory capacity exceeded.\n",
|
||||
"17940 Earth Engine memory capacity exceeded.\n",
|
||||
"17941 Earth Engine memory capacity exceeded.\n",
|
||||
"17945 Earth Engine memory capacity exceeded.\n",
|
||||
"17956 Earth Engine memory capacity exceeded.\n",
|
||||
"17958 Earth Engine memory capacity exceeded.\n",
|
||||
"17967 Earth Engine memory capacity exceeded.\n",
|
||||
"17971 Earth Engine memory capacity exceeded.\n",
|
||||
"17999 Earth Engine memory capacity exceeded.\n",
|
||||
"18012 Earth Engine memory capacity exceeded.\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import time \n",
|
||||
"cached_ids = get_cached_ids()\n",
|
||||
"\n",
|
||||
"def get_image_for_leak(i, cached_ids=cached_ids): \n",
|
||||
" leak = leaks_ATX.iloc[[i]]\n",
|
||||
" repo_date_ts = arrow.get(leak.REPO_Date.values[0]).timestamp\n",
|
||||
"# distance = resolution_min*(pixel_length/2.0+1)\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" # crappy way or recording that we tried this one\n",
|
||||
" leak_id = str(leak.OBJECTID.values[0])\n",
|
||||
" if leak_id in cached_ids:\n",
|
||||
" logger.info('Skipping cached download for leak id %s ',leak_id)\n",
|
||||
" return\n",
|
||||
" \n",
|
||||
" boundary = get_boundary(leak, distance=distance)\n",
|
||||
" \n",
|
||||
" # get image day before \n",
|
||||
" sentinel2_before = ee.ImageCollection('LANDSAT/LE7_L1T')\\\n",
|
||||
" .filterBounds(boundary)\\\n",
|
||||
" .filterDate((repo_date_ts-time_bin_delta)*1000,(repo_date_ts)*1000)\\\n",
|
||||
" .sort('system:time_start', opt_ascending=False) # first will be latest\n",
|
||||
" \n",
|
||||
" results = sentinel2_before.size().getInfo()\n",
|
||||
" if results<1:\n",
|
||||
" logger.info('Error no results for day before %s',leak_id)\n",
|
||||
" cached_ids = init_cache(leak_id) # so we know there where no results\n",
|
||||
" return\n",
|
||||
" \n",
|
||||
" # get image day after\n",
|
||||
" sentinel2_after = ee.ImageCollection('LANDSAT/LE7_L1T')\\\n",
|
||||
" .filterBounds(boundary)\\\n",
|
||||
" .filterDate((repo_date_ts)*1000,(repo_date_ts+time_bin_delta*6)*1000)\\\n",
|
||||
" .sort('system:time_start', opt_ascending=True) # first will be earliest\n",
|
||||
" \n",
|
||||
" results = sentinel2_after.size().getInfo()\n",
|
||||
" if results<1:\n",
|
||||
" logger.info('Error no results for day after, id %s',leak_id)\n",
|
||||
" cached_ids = init_cache(leak_id) # so we know there where no results\n",
|
||||
" return\n",
|
||||
" \n",
|
||||
" # download as save images \n",
|
||||
" logger.info('results for %s', leak_id)\n",
|
||||
" image = ee.Image(sentinel2_before.first()).clip(boundary)\n",
|
||||
" name=leak_id+'_before'\n",
|
||||
" path,files=download_image(\n",
|
||||
" image, \n",
|
||||
" scale=resolution_min, \n",
|
||||
" crs=crs_grid, \n",
|
||||
" name=name,\n",
|
||||
" cache_dir=cache_dir\n",
|
||||
" )\n",
|
||||
" # also save metadata so we can filter by date\n",
|
||||
" with open(path.joinpath('metadata.json'), 'w') as fo:\n",
|
||||
" metadata = dict(\n",
|
||||
" image=image.getInfo(),\n",
|
||||
" scale=resolution_min,\n",
|
||||
" crs=crs_grid,\n",
|
||||
" name=name,\n",
|
||||
" distance=distance,\n",
|
||||
" leak=json.loads(leak.to_json())\n",
|
||||
" )\n",
|
||||
" json.dump(metadata, fo)\n",
|
||||
"\n",
|
||||
" image = ee.Image(sentinel2_after.first()).clip(boundary)\n",
|
||||
" name=leak_id+'_after'\n",
|
||||
" path,files=download_image(\n",
|
||||
" image, \n",
|
||||
" scale=resolution_min, \n",
|
||||
" crs=crs_grid, \n",
|
||||
" name=name,\n",
|
||||
" cache_dir=cache_dir\n",
|
||||
" )\n",
|
||||
" with open(path.joinpath('metadata.json'), 'w') as fo:\n",
|
||||
" metadata = dict(\n",
|
||||
" image=image.getInfo(),\n",
|
||||
" scale=resolution_min,\n",
|
||||
" crs=crs_grid,\n",
|
||||
" name=name,\n",
|
||||
" distance=distance,\n",
|
||||
" leak=json.loads(leak.to_json())\n",
|
||||
" )\n",
|
||||
" json.dump(metadata, fo)\n",
|
||||
" \n",
|
||||
"for i in tqdm(range(len(leaks_ATX))):\n",
|
||||
" try:\n",
|
||||
" get_image_for_leak(i)\n",
|
||||
" except urllib.error.HTTPError as e:\n",
|
||||
" print(i,e)\n",
|
||||
" if e.code == 429:\n",
|
||||
" time.sleep(13);\n",
|
||||
" except Exception as e:\n",
|
||||
" print(i,e)\n",
|
||||
" ee.Initialize()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-18T04:36:19.085524Z",
|
||||
"start_time": "2017-03-18T12:36:19.082804+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# parsing tiffs"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-15T04:39:05.795321Z",
|
||||
"start_time": "2017-03-15T12:39:05.788009+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-20T06:54:54.295929Z",
|
||||
"start_time": "2017-03-20T13:49:03.076436+08:00"
|
||||
},
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"application/vnd.jupyter.widget-view+json": {
|
||||
"model_id": "b66a8a0c73ef43569f24d1934b59e24a"
|
||||
}
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Invalid metadata.json, deleted folder ../data/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2/ee_l7_AUTX-leaks_cache_v2/65498_before_3857_15.0, please rerun scraping cell to rescrape this image\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(40745, 0)"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# This loads it as X and y for machine learning, and also time and metadata so we can filter\n",
|
||||
"import shapely\n",
|
||||
"X = []\n",
|
||||
"y = []\n",
|
||||
"t = []\n",
|
||||
"m = []\n",
|
||||
"discarded=[]\n",
|
||||
"cdirs = [cdir for cdir in cache_dir.listdir() if ('_after_' in cdir) or ('_before_' in cdir)]\n",
|
||||
"for path in tqdm(cdirs):\n",
|
||||
" if not path.isdir(): continue\n",
|
||||
" files = [file.relpath(path) for file in path.listdir() if file.isfile() and file.endswith('.tif')]\n",
|
||||
" if files:\n",
|
||||
" # check metadata\n",
|
||||
" try:\n",
|
||||
" metadata = json.load(open(path.joinpath('metadata.json')))\n",
|
||||
" except (FileNotFoundError, ValueError) as e:\n",
|
||||
" path.move(path.dirname().dirname().joinpath('.deleteme-'+str(uuid.uuid4())))\n",
|
||||
" if '_after_' in path: # also delete the before path \n",
|
||||
" path_after = Path(path.replace('_after_','_before_'))\n",
|
||||
" if path_after.isdir():\n",
|
||||
" path_after.move(path.dirname().dirname().joinpath('.deleteme-'+str(uuid.uuid4())))\n",
|
||||
" logger.error('Invalid metadata.json, deleted folder %s, please rerun scraping cell to rescrape this image', path)\n",
|
||||
" continue\n",
|
||||
" \n",
|
||||
" # e.g. lets filter it so \"before\" image are only 1 day before\n",
|
||||
" if '_before_' in path.basename():\n",
|
||||
" yy = True\n",
|
||||
" else:\n",
|
||||
" yy = False\n",
|
||||
" \n",
|
||||
" # work out time gap too\n",
|
||||
" t1 = arrow.get(metadata['image']['properties']['system:time_end']/1000)\n",
|
||||
" t0 = arrow.get(metadata['leak']['features'][0]['properties']['REPO_Date'])\n",
|
||||
" td=t1-t0\n",
|
||||
" tt = td.total_seconds()\n",
|
||||
" \n",
|
||||
" # load data\n",
|
||||
" data = tifs2np(path,files,bands=bands)\n",
|
||||
" \n",
|
||||
" # check we don't have empty bands 1-13\n",
|
||||
" empty_bands = np.array([d.sum() for d in data])==0\n",
|
||||
" \n",
|
||||
" # lets check we didn't get the edge of an image\n",
|
||||
" bbox = np.array(metadata['image']['properties']['system:footprint']['coordinates'][0])\n",
|
||||
" loc = metadata['leak']['features'][0]['geometry']['coordinates']\n",
|
||||
" minx=bbox[:,0].min()\n",
|
||||
" maxx=bbox[:,0].max()\n",
|
||||
" miny=bbox[:,1].min()\n",
|
||||
" maxy=bbox[:,1].max()\n",
|
||||
" bbox_shp = shapely.geometry.box(\n",
|
||||
" minx=minx,\n",
|
||||
" maxx=maxx,\n",
|
||||
" miny=miny,\n",
|
||||
" maxy=maxy\n",
|
||||
" )\n",
|
||||
" loc_shp = shapely.geometry.Point(loc[0],loc[1])\n",
|
||||
" shapely.geometry.GeometryCollection([bbox_shp, loc_shp])\n",
|
||||
" try:\n",
|
||||
" assert loc_shp.intersects(bbox_shp), 'leak location should be inside image'\n",
|
||||
" assert bbox_shp.centroid.almost_equals(loc_shp, decimal=5), 'leak should be near center of image'\n",
|
||||
" assert (np.array([d.shape for d in data])==pixel_length).all(), 'image area should be the right amount of pixels'\n",
|
||||
" assert (maxx-minx)/(maxy-miny)<1.3, 'should be roughly square'\n",
|
||||
" assert (maxx-minx)/(maxy-miny)>0.7, 'should be roughly square'\n",
|
||||
" assert not empty_bands.all(), 'should not have all bands empty'\n",
|
||||
" except Exception as exc:\n",
|
||||
" print(path, exc)\n",
|
||||
"# raise(exc)\n",
|
||||
" discarded.append(path)\n",
|
||||
" else:\n",
|
||||
" X.append(data)\n",
|
||||
" y.append(yy)\n",
|
||||
" t.append(tt)\n",
|
||||
" m.append(metadata)\n",
|
||||
"\n",
|
||||
"len(X), len(discarded)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-20T05:48:19.175469Z",
|
||||
"start_time": "2017-03-20T13:48:07.994632+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-20T07:48:45.188727Z",
|
||||
"start_time": "2017-03-20T15:48:45.173290+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Path('../data/20170314-05-26-52_testing_earth_engine-l7-AUTX_v2/ee_l7_AUTX-leaks_cache_v2/70616_before_3857_15.0')"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"path"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"end_time": "2017-03-20T07:48:45.456224Z",
|
||||
"start_time": "2017-03-20T15:48:45.389047+08:00"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# shuffle\n",
|
||||
"from sklearn.utils import shuffle\n",
|
||||
"X,y,m= shuffle(X,y,m,random_state=1337)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"start_time": "2017-03-20T07:48:45.667Z"
|
||||
}
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"amount of each band [('B1', 0), ('B2', 0), ('B3', 0), ('B4', 0), ('B5', 0), ('B6_VCID_1', 0), ('B6_VCID_2', 0), ('B7', 0), ('B8', 0)]\n",
|
||||
"mean amount of bands 0.0\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# which bands do we have?\n",
|
||||
"a=np.array([x.sum(-1).sum(-1)==0 for x in X])\n",
|
||||
"print('amount of each band',list(zip(bands,a.sum(0))))\n",
|
||||
"print('mean amount of bands',a.sum(1).mean())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"start_time": "2017-03-20T07:48:46.082Z"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# save using hdf5 (so keras can easily load it) and json \n",
|
||||
"import h5py\n",
|
||||
"h5file = output_dir.joinpath('data.h5')\n",
|
||||
"with h5py.File(h5file, 'w') as h5f:\n",
|
||||
" h5f.create_dataset('X', data=X)\n",
|
||||
" h5f.create_dataset('y', data=y)\n",
|
||||
"\n",
|
||||
"json.dump(m,open(output_dir.joinpath('data_metadata.json'),'w'))\n",
|
||||
"\n",
|
||||
"with open(output_dir.joinpath('readme.md'),'w') as fo:\n",
|
||||
" fo. write(\"\"\"\n",
|
||||
"Files:\n",
|
||||
"- ee_S1_AUTX-leaks_cache- cached tiff files\n",
|
||||
"- script_metadata.json - information on scraping script\n",
|
||||
"- data.h5 contains X, y, and t.\n",
|
||||
" - X: tiff files for each band loaded into an array of shape (Leak, Bands, width, length)\n",
|
||||
" - y: True for before the leak, False for after\n",
|
||||
"- data_metadata: array of metadata for each leak in X. Each contain info on leak, image, and image search\n",
|
||||
" \n",
|
||||
"Loading: \n",
|
||||
"```py\n",
|
||||
"# load\n",
|
||||
"metadatas = json.load(open('data_metadata.json'))\n",
|
||||
"with h5py.File('data.h5','r') as h5f:\n",
|
||||
" X2 = h5f['X'][:]\n",
|
||||
" y2 = h5f['y'][:]\n",
|
||||
"y\n",
|
||||
"```\n",
|
||||
" \"\"\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"ExecuteTime": {
|
||||
"start_time": "2017-03-20T07:48:46.458Z"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# test load\n",
|
||||
"metadatas = json.load(open(output_dir.joinpath('data_metadata.json')))\n",
|
||||
"with h5py.File(output_dir.joinpath('data.h5'),'r') as h5f:\n",
|
||||
" X2 = h5f['X'][:]\n",
|
||||
" y2 = h5f['y'][:]\n",
|
||||
"X2.shape, y2, metadatas[0].keys()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (with sys packages)",
|
||||
"language": "python",
|
||||
"name": "py3syspck"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.4.2"
|
||||
},
|
||||
"toc": {
|
||||
"nav_menu": {
|
||||
"height": "96px",
|
||||
"width": "251px"
|
||||
},
|
||||
"navigate_menu": true,
|
||||
"number_sections": true,
|
||||
"sideBar": true,
|
||||
"threshold": 4,
|
||||
"toc_cell": false,
|
||||
"toc_section_display": "block",
|
||||
"toc_window_display": false
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 1
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
+1034
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,184 @@
|
||||
# Leak detection
|
||||
|
||||
*Michael Clark, Craig Baldachinno*
|
||||
|
||||
This project explores the idea that we can find leaks using satellite images and machine learning. The primary hypothesis is this: given the time and place of leak repairs and satellite images of 10-15m2 resolution we predict leaks better than a baseline defined by random guessing. Our results give and f1-score of around 0.6 compared to a dummy result of 0.5. This small result supports our hypothesis but the predictive capability is too small to be used for pipe maintenance. Our conclusion is that we can make small and uncertain predictions but with higher resolution data predictive capability may improve to useful amounts.
|
||||
|
||||
# Introduction
|
||||
|
||||
Water leakage in Western Australia account for more than 10 billion litres of water and costs over a billion dollars each year. The Water Corporation can only manually survey 10-12 percent of the pipe network each year for non-visible leaks. [[0]](https://www.watercorporation.com.au/-/media/files/about-us/our-performance/annual-report-2015/water-corporation-annual-report-2015.pdf). A solution that can detect leaks without additional on-site inspections would be valuable for water utilities around the world.
|
||||
|
||||
One solution is to use satellite images. Satellite derived vegetation indices can be used for the detection of water pipeline leakages in semiarid areas. [[1]](http://spie.org/Publications/Proceedings/Paper/10.1117/12.2028241). Spectral signiture profiles can also be used to detect leaks but the signitures can be complicated [[2]](http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1757291).
|
||||
|
||||
Deep learning has been successful at detecting objects from satellite images even when the images can give complicated response profiles [[3]](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.232.1679&rep=rep1&type=pdf). This project attempts to find leaks using satellite images - utilizing machine learning to overcome difficulties with compelex data.
|
||||
|
||||
# Data
|
||||
|
||||
We have a leak datasets, obtained by a freedom of information request to the very helpful people at Austin Water. Quality is hard to verify since little context is available. For example some entries could be meter repairs, false call-outs, typos, duplicates, or multi-stage repairs.
|
||||
|
||||
|
||||
These have been combined with satellite image from Google Earth Engine to make machine learning datasets. For more see [the scraping notebooks](notebooks/scraped_satellite_images) and the readme files for each dataset. An example of this data is below:
|
||||
|
||||
|
||||
|
||||
Satellites only pass over once a week at best so not every leak has a corresponding image.
|
||||
|
||||
|
||||
| dataset | Satellite | Leak dataset | Resolution ** |
|
||||
|-----------------------------|--------------|----------------|-----------------|
|
||||
| s2-AUTX_v6_COPERNICUS-S2 | [Sentinel 2](https://explorer.earthengine.google.com/#detail/COPERNICUS%2FS2) | ATX | >10m2 |
|
||||
| l8-AUTX_v2_LANDSAT-LC8_L1T | [Landsat 8](https://explorer.earthengine.google.com/#detail/LANDSAT%2FLC8_L1T) | ATX | >15m2 |
|
||||
| l7-AUTX_v2_LANDSAT-LE7_L1T | [Landsat 7](https://explorer.earthengine.google.com/#detail/LANDSAT%2FLE7_L1T) | ATX | >15m2 |
|
||||
| s1-AURX_COPERNICUS-S1_GRD | [Sentinel 1](https://explorer.earthengine.google.com/#detail/COPERNICUS%2FS1_GRD) | ATX | >10m2 |
|
||||
|
||||
Notes:
|
||||
|
||||
\* This data is combined from other rows in the table
|
||||
|
||||
** This gives source resolution. All data was interpolated
|
||||
|
||||
What data filters yield the best results? To answer this question we used hyperparameter optimisation to find the data filters that have the most predictive power for each dataset. For more see the [hyperopt notebook](notebooks/3a_hyperopt.ipynb) and the data tables in the appendix.
|
||||
|
||||
# Models
|
||||
|
||||
Modelling and prediction is done in the [3_model notebook](notebooks/3_model.ipynb) where outputs include a map of false positives, predictions, and metrics. We tried many models and archetectures including U-Net (citation), pretrained VGG-16 (citation), decision-tree regressor, multi-layer perceptrons. No model gave any significant improvement over random forest methods so this was our primary benchmark.
|
||||
|
||||
# Results
|
||||
|
||||
Using a random forest model we tested 500 models on each of our satellite datasets. Dummy scores were generated by scikit's dummy classifiers with a uniform strategy. Since we are using balanced data we can use f1-scores as our primary benchmark. If max score doesn't exceed the dummy score then it's a dataset with lots of misleading features. The table below shows the results, for more see the [2_qc_scraped_data notebook](notebooks/2_qc_scraped_data.ipynb).
|
||||
|
||||
| | support | dummy_f1_leak | f1_leak_max | f1_leak_mean |
|
||||
|:----------------------------|----------:|----------------:|--------------:|---------------:|
|
||||
| s2-AUTX_COPERNICUS-S2 | 224 | 0.57 | 0.66 | 0.521 |
|
||||
| l7-AUTX_LANDSAT-LE7_L1T | 1082 | 0.54 | 0.66 | 0.426 |
|
||||
| l8-AUTX_LANDSAT-LC8_L1T | 434 | 0.55 | 0.69 | 0.55 |
|
||||
| s1-AUTX_COPERNICUS-S1_GRD | 281 | 0.58 | 0.68 | 0.511 |
|
||||
|
||||
Our results show that random forest can improve on dummy benchmarks. However the difference is small, and the f1-scores are low, implying to many false positives to be useful to water utility repair teams.
|
||||
|
||||
# Conclusion
|
||||
|
||||
The results support our hypothesis, that machine learning can use satellite images (of 10m2 resolution) to predict water leaks. However the f1-scores indicate predictions that are to unreliable to be used in the field.
|
||||
|
||||
The largest limitation in this experiment was the resolution of satellite images. We expect better results with higher resolution images, unfortunately only a limited amout of high resolution imagery is available in 2017. Using the entire digital globe catalogue would yield around 20 matches with our leaks datasets, and the same goes for the high resolution US National Agriculture Imagery Program (NAIP).
|
||||
|
||||
|
||||
# Appendix
|
||||
|
||||
## Files
|
||||
|
||||
The main files are
|
||||
- `./notebooks/0_process_data` where the leak datasets are converted
|
||||
- `./notebooks/scraped_satellite_images` where images are downloaded for each leak
|
||||
- `notebooks/2_qc_scraped_data.ipynb` where test each dataset for consistency and predictive power
|
||||
- `notebooks/2a_hyperopt.ipynb` where I try varied filters to see which has the most predictive power
|
||||
- `./notebooks/2_model.ipynb` which is a jupyter-notebook script for python 3.
|
||||
|
||||
## Data filters
|
||||
|
||||
We have many options for filtering the data including balancing datasets, jittering value, augumenting data, and filtering out cloudy images. The most important parameter is "timespan_before", here we decide how long we will allow between a satellite image and a leak repair. If we choose to long we will have a image of a pipe before it started leaking, and if we choose to short a time we will reduce our dataset.
|
||||
|
||||
To solve this I ran hyperparamater optimisation where I tried all variations using a random forest model. The results are below. The results for Landsat-7 show that the best 50 models used a timespan_before of 1.13 days, before which images are discarded.
|
||||
|
||||
### best Landsat-7 filters
|
||||
|
||||
top 50 scores,
|
||||
|
||||
- by metric="result_metrics_report_f1-score_leak"
|
||||
- for key="../output/hyperopt/data_filters_derived_l7_ATX.json"
|
||||
- number of trials=564
|
||||
|
||||
| | corr | mean | std | units |
|
||||
|:------------------------------------------------|--------:|--------:|--------:|:--------|
|
||||
| misc_vals_balanced_classes | 0.118 | 0.98 | 0.141 | bool |
|
||||
| misc_vals_channel_shift_range | -0.137 | 0.083 | 0.06 | frac |
|
||||
| misc_vals_height_shift_range | -0.138 | 0.165 | 0.094 | frac |
|
||||
| misc_vals_max_cloud_cover | -0.128 | 0.109 | 0.085 | frac |
|
||||
| misc_vals_normalized | nan | 0 | 0 | bool |
|
||||
| misc_vals_rescale | -0.011 | 0.33 | 0.089 | frac |
|
||||
| misc_vals_rotation_range | -0.041 | 9.903 | 7.671 | deg |
|
||||
| misc_vals_timespan_before | -0.24 | 1.134 | 1.181 | days |
|
||||
| misc_vals_width_shift_range | 0.027 | 0.271 | 0.091 | frac |
|
||||
| result_metrics_report_f1-score_leak | 1 | 0.679 | 0.01 | |
|
||||
| result_dummy_metrics_report_dummy_f1-score_leak | 0.218 | 0.508 | 0.039 | |
|
||||
| result_dummy_metrics_report_dummy_support_leak | -0.259 | 217.64 | 207.412 | int |
|
||||
|
||||
|
||||
Here's how I interpret this table.
|
||||
|
||||
- balanced_classes: For the l7 datasets should we used balanced classes?
|
||||
- We have a decent positive correlation so it's important.
|
||||
- We had a mean of 0.98, but it's bool, so we should set it to True.
|
||||
|
||||
- channel_shift_range: How about channel shift range?
|
||||
- This has a decent negative correlation so it's important and having it too high gives a low f1 score.
|
||||
- The mean is 0.08+-0.06 so lets set it to 0.8.
|
||||
|
||||
### Best Landsat-8 filters
|
||||
|
||||
top 50 scores,
|
||||
|
||||
- by metric="result_metrics_report_f1-score_leak"
|
||||
- for key="../output/hyperopt/data_filters_derived_l8_joined_ATX.json"
|
||||
- number of trials=564
|
||||
|
||||
| | corr | mean | std | units |
|
||||
|:------------------------------------------------|-------:|--------:|--------:|:--------|
|
||||
| misc_vals_balanced_classes | 0.033 | 0.6 | 0.495 | bool |
|
||||
| misc_vals_channel_shift_range | 0.055 | 0.107 | 0.051 | frac |
|
||||
| misc_vals_height_shift_range | -0.04 | 0.22 | 0.147 | frac |
|
||||
| misc_vals_max_cloud_cover | -0.116 | 0.25 | 0.257 | frac |
|
||||
| misc_vals_normalized | -0.784 | 0.42 | 0.499 | bool |
|
||||
| misc_vals_rescale | 0.292 | 0.15 | 0.107 | frac |
|
||||
| misc_vals_rotation_range | -0.072 | 21.022 | 9.381 | deg |
|
||||
| misc_vals_timespan_before | -0.243 | 5.069 | 1.431 | days |
|
||||
| misc_vals_width_shift_range | 0.082 | 0.163 | 0.097 | frac |
|
||||
| result_metrics_report_f1-score_leak | 1 | 0.633 | 0.027 | |
|
||||
| result_dummy_metrics_report_dummy_f1-score_leak | 0.006 | 0.495 | 0.026 | |
|
||||
| result_dummy_metrics_report_dummy_support_leak | -0.235 | 411.2 | 173.273 | int |
|
||||
|
||||
### Best Sentinel-2 filters
|
||||
|
||||
top 50 scores,
|
||||
|
||||
- by metric="result_metrics_report_f1-score_leak"
|
||||
- for key="../output/hyperopt/derived_s2_joined_ATX.json"
|
||||
- number of trials=1000
|
||||
|
||||
| | mean | std | corr | units |
|
||||
|:------------------------------------------------|--------:|-------:|-------:|:--------|
|
||||
| misc_vals_balanced_classes | 0.98 | 0.141 | 0.027 | bool |
|
||||
| misc_vals_channel_shift_range | 0.129 | 0.071 | -0.087 | frac |
|
||||
| misc_vals_height_shift_range | 0.35 | 0.134 | 0.008 | frac |
|
||||
| misc_vals_max_cloud_cover | 0.171 | 0.168 | -0.222 | frac |
|
||||
| misc_vals_normalized | 0.28 | 0.454 | 0.2 | bool |
|
||||
| misc_vals_rescale | 0.142 | 0.106 | -0.063 | frac |
|
||||
| misc_vals_rotation_range | 21.399 | 13.651 | -0.092 | deg |
|
||||
| misc_vals_timespan_before | 1.984 | 1.769 | -0.238 | days |
|
||||
| misc_vals_width_shift_range | 0.418 | 0.099 | -0.302 | frac |
|
||||
| result_metrics_report_f1-score_leak | 0.684 | 0.021 | 1 | |
|
||||
| result_dummy_metrics_report_dummy_f1-score_leak | 0.505 | 0.061 | 0.105 | |
|
||||
| result_dummy_metrics_report_dummy_support_leak | 113.68 | 88.334 | -0.35 | int |
|
||||
|
||||
### Best Sentinel-1 filters
|
||||
|
||||
top 50 scores,
|
||||
|
||||
- by metric="result_metrics_report_f1-score_leak"
|
||||
- for key="../output/hyperopt/data_filters_s1-AUTX_COPERNICUS-S1_GRD.json"
|
||||
- number of trials=1000
|
||||
|
||||
| | mean | std | corr | units |
|
||||
|:------------------------------------------------|--------:|--------:|--------:|:--------|
|
||||
| misc_vals_balanced_classes | 1 | 0 | nan | bool |
|
||||
| misc_vals_channel_shift_range | 0.169 | 0.06 | -0.063 | frac |
|
||||
| misc_vals_height_shift_range | 0.167 | 0.115 | 0.146 | frac |
|
||||
| misc_vals_max_cloud_cover | 0.485 | 0.286 | -0.121 | frac |
|
||||
| misc_vals_normalized | 0.06 | 0.24 | -0.028 | bool |
|
||||
| misc_vals_rescale | 0.266 | 0.125 | -0.105 | frac |
|
||||
| misc_vals_rotation_range | 19.203 | 9.334 | 0.232 | deg |
|
||||
| misc_vals_timespan_before | 2.099 | 1.21 | -0.285 | days |
|
||||
| misc_vals_width_shift_range | 0.261 | 0.111 | 0.075 | frac |
|
||||
| result_metrics_report_f1-score_leak | 0.668 | 0.014 | 1 | |
|
||||
| result_dummy_metrics_report_dummy_f1-score_leak | 0.508 | 0.04 | 0.093 | |
|
||||
| result_dummy_metrics_report_dummy_support_leak | 165.88 | 100.096 | -0.237 | int |
|
||||
Reference in New Issue
Block a user