From eecdb485fc28c202f8fa7db85323ce3f2881c24d Mon Sep 17 00:00:00 2001 From: Kilian Vos Date: Tue, 3 Apr 2018 12:03:38 +1000 Subject: [PATCH] new implementation --- download_images.py | 118 ++++++++++++++++++++++++++++++++++++++ read_images.py | 140 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 258 insertions(+) create mode 100644 download_images.py create mode 100644 read_images.py diff --git a/download_images.py b/download_images.py new file mode 100644 index 0000000..59279b6 --- /dev/null +++ b/download_images.py @@ -0,0 +1,118 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Mar 27 17:12:35 2018 + +@author: Kilian +""" + +# Initial settings +import os +import numpy as np +import matplotlib.pyplot as plt +import pdb +import ee + + +# other modules +from osgeo import gdal, ogr, osr +from urllib.request import urlretrieve +import zipfile +from datetime import datetime +import pytz +import pickle + + + +# 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.measure as measure + + +# import own modules +from functions.utils import * + +np.seterr(all='ignore') # raise/ignore divisions by 0 and nans +ee.Initialize() + +def download_tif(image, polygon, bandsId, filepath): + """downloads tif image (region and bands) from the ee server and stores it in a temp file""" + url = ee.data.makeDownloadUrl(ee.data.getDownloadId({ + 'image': image.serialize(), + 'region': polygon, + 'bands': bandsId, + 'filePerBand': 'false', + 'name': 'data', + })) + local_zip, headers = urlretrieve(url) + with zipfile.ZipFile(local_zip) as local_zipfile: + return local_zipfile.extract('data.tif', filepath) + +# select collection +input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA') +# location (Narrabeen-Collaroy beach) +rect_narra = [[[151.301454, -33.700754], + [151.311453, -33.702075], + [151.307237, -33.739761], + [151.294220, -33.736329], + [151.301454, -33.700754]]]; + +# dates +start_date = '2016-01-01' +end_date = '2016-12-31' +# filter by location +flt_col = input_col.filterBounds(ee.Geometry.Polygon(rect_narra)).filterDate(start_date, end_date) + +n_img = flt_col.size().getInfo() +print('Number of images covering Narrabeen:', n_img) +im_all = flt_col.getInfo().get('features') + +satname = 'L8' +sitename = 'NARRA' +suffix = '.tif' +filepath = os.path.join(os.getcwd(), 'data', satname, sitename) +filepath_pan = os.path.join(filepath, 'pan') +filepath_ms = os.path.join(filepath, 'ms') + +all_names_pan = [] +all_names_ms = [] +timestamps = [] +# loop through all images +for i in range(n_img): + # find each image in ee database + im = ee.Image(im_all[i].get('id')) + im_dic = im.getInfo() + im_bands = im_dic.get('bands') + im_date = im_dic['properties']['DATE_ACQUIRED'] + t = im_dic['properties']['system:time_start'] + im_timestamp = datetime.fromtimestamp(t/1000, tz=pytz.utc) + timestamps.append(im_timestamp) + im_epsg = int(im_dic['bands'][0]['crs'][5:]) + + # delete dimensions key from dictionnary, otherwise the entire image is extracted + for j in range(len(im_bands)): del im_bands[j]['dimensions'] + pan_band = [im_bands[7]] + ms_bands = [im_bands[1], im_bands[2], im_bands[3], im_bands[4], im_bands[5], im_bands[11]] + + filename_pan = satname + '_' + sitename + '_' + im_date + '_pan' + suffix + filename_ms = satname + '_' + sitename + '_' + im_date + '_ms' + suffix + + print(i) + if any(filename_pan in _ for _ in all_names_pan): + filename_pan = satname + '_' + sitename + '_' + im_date + '_pan' + '_r' + suffix + filename_ms = satname + '_' + sitename + '_' + im_date + '_ms' + '_r' + suffix + all_names_pan.append(filename_pan) + +# local_data_pan = download_tif(im, rect_narra, pan_band, filepath_pan) +# os.rename(local_data_pan, os.path.join(filepath_pan, filename_pan)) +# local_data_ms = download_tif(im, rect_narra, ms_bands, filepath_ms) +# os.rename(local_data_ms, os.path.join(filepath_ms, filename_ms)) + + +with open(os.path.join(filepath, sitename + '_timestamps' + '.pkl'), 'wb') as f: + pickle.dump(timestamps, f) +with open(os.path.join(filepath, sitename + '_epsgcode' + '.pkl'), 'wb') as f: + pickle.dump(im_epsg, f) + \ No newline at end of file diff --git a/read_images.py b/read_images.py new file mode 100644 index 0000000..4b08ff9 --- /dev/null +++ b/read_images.py @@ -0,0 +1,140 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Mar 27 17:12:35 2018 + +@author: Kilian +""" + +# Initial settings +import os +import numpy as np +import matplotlib.pyplot as plt +import ee +import pdb + +# other modules +from osgeo import gdal, ogr, osr +import pickle +import matplotlib.cm as cm +from pylab import ginput + +# 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.measure as measure + + +# import own modules +import functions.utils as utils +import functions.sds as sds + +np.seterr(all='ignore') # raise/ignore divisions by 0 and nans +ee.Initialize() + +# initial settings +cloud_thresh = 0.5 # threshold for cloud cover +plot_bool = False # 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 +output_epsg = 28356 # GDA94 / MGA Zone 56 + +satname = 'L8' +sitename = 'NARRA' +filepath = os.path.join(os.getcwd(), 'data', satname, sitename) +with open(os.path.join(filepath, sitename + '_timestamps' + '.pkl'), 'rb') as f: + timestamps = pickle.load(f) +timestamps_sorted = sorted(timestamps) +with open(os.path.join(filepath, sitename + '_epsgcode' + '.pkl'), 'rb') as f: + input_epsg = pickle.load(f) + + +file_path_pan = os.path.join(os.getcwd(), 'data', 'L8', 'NARRA', 'pan') +file_path_ms = os.path.join(os.getcwd(), 'data', 'L8', 'NARRA', 'ms') +file_names_pan = os.listdir(file_path_pan) +file_names_ms = os.listdir(file_path_ms) +N = len(file_names_pan) +idx_high_cloud = [] +t = [] +shorelines = [] + +for i in range(N): + # read pan image + fn_pan = os.path.join(file_path_pan, file_names_pan[i]) + data = gdal.Open(fn_pan, gdal.GA_ReadOnly) + georef = np.array(data.GetGeoTransform()) + bands = [data.GetRasterBand(i + 1).ReadAsArray() for i in range(data.RasterCount)] + im_pan = np.stack(bands, 2)[:,:,0] + # read ms image + fn_ms = os.path.join(file_path_ms, file_names_ms[i]) + data = gdal.Open(fn_ms, gdal.GA_ReadOnly) + bands = [data.GetRasterBand(i + 1).ReadAsArray() for i in range(data.RasterCount)] + im_ms = np.stack(bands, 2) + # cloud mask + im_qa = im_ms[:,:,5] + cloud_mask = sds.create_cloud_mask(im_qa) + cloud_mask = transform.resize(cloud_mask, (im_pan.shape[0], im_pan.shape[1]), + order=0, preserve_range=True, + mode='constant').astype('bool_') + # resize the image using bilinear interpolation (order 1) + im_ms = transform.resize(im_ms,(im_pan.shape[0], im_pan.shape[1]), + order=1, preserve_range=True, mode='constant') + # check if -inf or nan values and add to cloud mask + im_inf = np.isin(im_ms[:,:,0], -np.inf) + im_nan = np.isnan(im_ms[:,:,0]) + cloud_mask = np.logical_or(np.logical_or(cloud_mask, im_inf), im_nan) + cloud_content = sum(sum(cloud_mask.astype(int)))/(cloud_mask.shape[0]*cloud_mask.shape[1]) + if cloud_content > cloud_thresh: + print('skipped ' + str(i)) + idx_high_cloud.append(i) + continue + # rescale intensities + im_ms = sds.rescale_image_intensity(im_ms, cloud_mask, prob_high, plot_bool) + im_pan = sds.rescale_image_intensity(im_pan, cloud_mask, prob_high, plot_bool) + # pansharpen rgb image + im_ms_ps = sds.pansharpen(im_ms[:,:,[0,1,2]], im_pan, cloud_mask, 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 = sds.nd_index(im_ms_ps[:,:,3], im_ms_ps[:,:,1], cloud_mask, plot_bool) + # detect edges + wl_pix = sds.find_wl_contours(im_ndwi, cloud_mask, min_contour_points, plot_bool) + # convert from pixels to world coordinates + wl_coords = sds.convert_pix2world(wl_pix, georef) + # convert to output epsg spatial reference + wl = sds.convert_epsg(wl_coords, input_epsg, output_epsg) + +# plt.figure() +# plt.imshow(im_ms_ps[:,:,[2,1,0]]) +# for i,contour in enumerate(wl_pix): plt.plot(contour[:, 1], contour[:, 0], linewidth=2) +# plt.axis('image') +# plt.title(file_names_pan[i]) +# plt.show() + + plt.figure() + centroids = [] + cmap = cm.get_cmap('jet') + for j,contour in enumerate(wl): + colours = cmap(np.linspace(0, 1, num=len(wl))) + centroids.append([np.mean(contour[:, 0]),np.mean(contour[:, 1])]) + plt.plot(contour[:, 0], contour[:, 1], linewidth=2, color=colours[j,:]) + plt.plot(np.mean(contour[:, 0]), np.mean(contour[:, 1]), 'o', color=colours[j,:]) + plt.axis('equal') + plt.title(file_names_pan[i]) + plt.draw() + pt_in = np.array(ginput(1)) + dist_centroid = [np.linalg.norm(_ - pt_in) for _ in centroids] + shorelines.append(wl[np.argmin(dist_centroid)]) + t.append(timestamps_sorted[i]) + +#plt.figure() +#plt.axis('equal') +#for j in range(len(shorelines)): +# plt.plot(shorelines[j][:,0], shorelines[j][:,1]) +#plt.draw() + +output = {'t':t, 'shorelines':shorelines} + +with open(os.path.join(filepath, sitename + '_output' + '.pkl'), 'wb') as f: + pickle.dump(output, f) \ No newline at end of file