{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport sklearn\nimport glob,os\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\nsmjpegs = [f for f in glob.glob(\"../input/draper-satellite-image-chronology/train_sm/train_sm/*.jpeg\")]\nprint(smjpegs[:9])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"set175 = [smj for smj in smjpegs if \"set175\" in smj]\nprint(set175)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first=plt.imread(\"../input/draper-satellite-image-chronology/train_sm/train_sm/set175_1.jpeg\")\ndims=np.shape(first)\nprint(dims)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.min(first),np.max(first)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pixel_matrix=np.reshape(first,(dims[0]* dims[1],dims[2]))\nprint(np.shape(pixel_matrix))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plt.scatter(pixel_matrix[:,0], pixel_matrix[:,1])\n_ = plt.hist2d(pixel_matrix[:,1], pixel_matrix[:,2], bins=(50,50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fifth=plt.imread(\"../input/draper-satellite-image-chronology/train_sm/train_sm/set175_1.jpeg\")\ndims=np.shape(fifth)\npixel_matrix5=np.reshape(fifth,(dims[0]* dims[1],dims[2]))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = plt.hist2d(pixel_matrix5[:,1], pixel_matrix5[:,2], bins=(50,50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = plt.hist2d(pixel_matrix[:,2], pixel_matrix5[:,2], bins=(50,50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(first);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(first[:,:,2] - fifth[:,:,1]);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"second=plt.imread(\"../input/draper-satellite-image-chronology/train_sm/train_sm/set137_4.jpeg\")\n\nplt.imshow(second);\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"second = plt.imread('../input/draper-satellite-image-chronology/train_sm/train_sm/set175_2.jpeg')\nplt.imshow(first[:,:,2] - second[:,:,2]);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(second);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import cluster\n\n\nkmeans=cluster.KMeans(5)\nclustered=kmeans.fit_predict(pixel_matrix)\n\n\ndims=np.shape(first)\nclustered_img=np.reshape(clustered,(dims[0],dims[1]))\nplt.imshow(clustered_img)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(first);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ind0, ind1, ind2, ind3 = [np.where(clustered == x)[0] for x in [0, 1, 2, 3]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from mpl_toolkits.mplot3d import Axes3D\n\nfig=plt.figure()\nax=fig.add_subplot(111,projection=\"3d\")\nplot_vals = [('r', 'o', ind0),\n             ('b', '^', ind1),\n             ('g', '8', ind2),\n             ('m', '*', ind3)]\n\nfor c, m, ind in plot_vals:\n    xs = pixel_matrix[ind, 0]\n    ys = pixel_matrix[ind, 1]\n    zs = pixel_matrix[ind, 2]\n    ax.scatter(xs, ys, zs, c=c, marker=m)\n\nax.set_xlabel('Blue channel')\nax.set_ylabel('green channel')\nax.set_zlabel('Red channel')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import seaborn as sns\nsns.distplot(pixel_matrix[:,0], bins=12)\nsns.distplot(pixel_matrix[:,1], bins=12)\nsns.distplot(pixel_matrix[:,2], bins=12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"set79=[smj for smj in smjpegs if \"set79\" in smj]\nprint(set79)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img79_1, img79_2, img79_3, img79_4, img79_5 = \\\n  [plt.imread(\"../input/draper-satellite-image-chronology/train_sm/train_sm/set79_\" + str(n) + \".jpeg\") for n in range(1, 6)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_list = (img79_1, img79_2, img79_3, img79_4, img79_5)\n\nplt.figure(figsize=(8,10))\nplt.imshow(img_list[0])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class MSImage():\n    \"\"\"Lightweight wrapper for handling image to matrix transforms. No setters,\n    main point of class is to remember image dimensions despite transforms.\"\"\"\n    \n    def __init__(self, img):\n        \"\"\"Assume color channel interleave that holds true for this set.\"\"\"\n        self.img = img\n        self.dims = np.shape(img)\n        self.mat = np.reshape(img, (self.dims[0] * self.dims[1], self.dims[2]))\n\n    @property\n    def matrix(self):\n        return self.mat\n        \n    @property\n    def image(self):\n        return self.img\n    \n    def to_flat_img(self, derived):\n        \"\"\"\"Use dims property to reshape a derived matrix back into image form when\n        derived image would only have one band.\"\"\"\n        return np.reshape(derived, (self.dims[0], self.dims[1]))\n    \n    def to_matched_img(self, derived):\n        \"\"\"\"Use dims property to reshape a derived matrix back into image form.\"\"\"\n        return np.reshape(derived, (self.dims[0], self.dims[1], self.dims[2]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"msi79_1 = MSImage(img79_1)\nprint(np.shape(msi79_1.matrix))\nprint(np.shape(msi79_1.img))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def bnormalize(mat):\n    \"\"\"much faster brightness normalization, since it's all vectorized\"\"\"\n    bnorm = np.zeros_like(mat, dtype=np.float32)\n    maxes = np.max(mat, axis=1)\n    bnorm = mat / np.vstack((maxes, maxes, maxes)).T\n    return bnorm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"bnorm = bnormalize(msi79_1.matrix)\nbnorm_img = msi79_1.to_matched_img(bnorm)\nplt.figure(figsize=(8,10))\nplt.imshow(bnorm_img)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"msi79_2 = MSImage(img79_2)\nbnorm79_2 = bnormalize(msi79_2.matrix)\nbnorm79_2_img = msi79_2.to_matched_img(bnorm79_2)\nplt.figure(figsize=(8,10))\nplt.imshow(bnorm79_2_img)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"msinorm79_1 = MSImage(bnorm_img)\nmsinorm79_2 = MSImage(bnorm79_2_img)\n\n_ = plt.hist2d(msinorm79_1.matrix[:,2], msinorm79_2.matrix[:,2], bins=(50,50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = plt.hist2d(msinorm79_1.matrix[:,1], msinorm79_2.matrix[:,1], bins=(50,50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = plt.hist2d(msinorm79_1.matrix[:,0], msinorm79_2.matrix[:,0], bins=(50,50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8,10))\nplt.imshow(img79_1)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.max(img79_1[:,:,0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(img79_1[:,:,0] > 230)\nplt.subplot(122)\nplt.imshow(img79_1)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(img79_2[:,:,0] > 230)\nplt.subplot(122)\nplt.imshow(img79_2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.min(bnorm79_2_img[:,:,0]))\nprint(np.max(bnorm79_2_img[:,:,0]))\nprint(np.mean(bnorm79_2_img[:,:,0]))\nprint(np.std(bnorm79_2_img[:,:,0]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(bnorm79_2_img[:,:,0] > 0.98)\nplt.subplot(122)\nplt.imshow(img79_2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(bnorm_img[:,:,0] > 0.98)\nplt.subplot(122)\nplt.imshow(img79_1)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow((bnorm79_2_img[:,:,0] > 0.9999) & \\\n           (bnorm79_2_img[:,:,1] < 0.9999) & \\\n           (bnorm79_2_img[:,:,2] < 0.9999))\nplt.subplot(122)\nplt.imshow(img79_2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(bnorm_img[:,:,0] > 0.995)\nplt.subplot(122)\nplt.imshow(img79_1)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,6))\nplt.subplot(121)\nplt.plot(bnorm_img[2000, 1000, :])\nplt.subplot(122)\nplt.plot(img79_1[2000, 1000, :])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import spatial\npixel=msi79_1.matrix[2000*1000, :]\nnp.shape(pixel)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def spectral_angle_mapper(pixel):\n    return lambda p2: spatial.distance.cosine(pixel, p2)\n\nmatch_pixel = np.apply_along_axis(spectral_angle_mapper(pixel),\n1, msi79_1.matrix)\n\nplt.figure(figsize=(10,6))\nplt.imshow(msi79_1.to_flat_img(match_pixel < 0.0000001))\n\ndef summary(mat):\n    print(\"Max: \", np.max(mat),\n          \"Min: \", np.min(mat),\n          \"Std: \", np.std(mat),\n          \"Mean: \", np.mean(mat))\n\nsummary(match_pixel)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"set144=[MSImage(plt.imread(smj)) for smj in smjpegs if \"set144\" in smj]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(set144[0].image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import skimage\nfrom skimage.feature import greycomatrix, greycoprops\nfrom skimage.filters import sobel","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig=plt.figure(figsize=(8,8))\nplt.imshow(sobel(set144[0].image[:750,:750,2]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage.filters import sobel_h\n\n# can also apply sobel only across one direction.\nfig = plt.figure(figsize=(8, 8))\nplt.imshow(sobel_h(set144[0].image[:750,:750,2]), cmap='BuGn')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.decomposition import PCA\n\npca = PCA(3)\npca.fit(set144[0].matrix)\nset144_0_pca = pca.transform(set144[0].matrix)\nset144_0_pca_img = set144[0].to_matched_img(set144_0_pca)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig=plt.figure(figsize=(8,8))\nplt.imshow(set144_0_pca_img[:,:,0],cmap=\"BuGn\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig=plt.figure()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}