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

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/yoba7/Geomatics/blob/master/clustering_sur_orthophotos.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="vxhar62Ss2iS" colab_type="text"
# # Récupérer orthophotos avec owslib
#
# [owslib](https://geopython.github.io/OWSLib/) est une librairie qui permet de récupérer des informations d'un serveur WMS, WFS, ...

# + id="QsVmKvtI0Smc" colab_type="code" colab={}
# !pip install owslib opencv-python

# + [markdown] id="Um3Pbq0ptYTq" colab_type="text"
# ## Import de la librairie

# + id="kCMYsTH80RUP" colab_type="code" colab={}
from owslib.wms import WebMapService

# + [markdown] id="vAVz6oggt1MF" colab_type="text"
# ## Connexion avec un serveur
#
# Ici, je me connecte sur un serveur d'images du Service Public de Wallonie et je vérifie quels sont les formats supportés.

# + id="KCtSMy2e0RUS" colab_type="code" colab={}
#wms = WebMapService('http://geoservices.wallonie.be/arcgis/services/IMAGERIE/ORTHO_2018/MapServer/WMSServer')
#wms = WebMapService('http://geoservices.wallonie.be/arcgis/services/IMAGERIE/ORTHO_1971/MapServer/WMSServer?request=GetCapabilities&service=WMS')
#wms = WebMapService('http://geoservices.wallonie.be/arcgis/services/IMAGERIE/ORTHO_2009_2010/MapServer/WMSServer')
wms = WebMapService('http://geoservices.wallonie.be/arcgis/services/IMAGERIE/ORTHO_2012_2013/MapServer/WMSServer')

# + id="E098jRYG0RUU" colab_type="code" outputId="8edabcc2-2006-40ec-ae79-1de7aa35af1a" colab={"base_uri": "https://localhost:8080/", "height": 34}
list(wms.contents)

# + id="iLZGzoQx0RUZ" colab_type="code" outputId="ca9ac002-eed3-474e-dc1d-62fdc2cc495f" colab={"base_uri": "https://localhost:8080/", "height": 34}
wms['0'].title

# + id="ONmRa6g2kINs" colab_type="code" outputId="aad69d97-3076-409f-9309-79b27b7d6feb" colab={"base_uri": "https://localhost:8080/", "height": 173}
wms.getOperationByName('GetMap').formatOptions

# + id="B--baxsU0RUd" colab_type="code" outputId="217499cb-21ab-4d4a-ab62-de98b00e342f" colab={"base_uri": "https://localhost:8080/", "height": 72}
wms['0'].styles

# + [markdown] id="0l_gshQPt3kn" colab_type="text"
# ## Getmap
#
# Récupérer une orthophoto du Châteur de La Hulpe au format JPG

# + id="2H7tnAQT0RUh" colab_type="code" colab={}
# Chateau de La Hulpe
lowerLeftCorner=[156350,158090]
upperRightCorner=[x+300 for x in lowerLeftCorner]
myBBox=(lowerLeftCorner[0],lowerLeftCorner[1],upperRightCorner[0],upperRightCorner[1])

# + id="IaKRTt380RUj" colab_type="code" colab={}
img = wms.getmap(layers=['0'],
                 srs='EPSG:31370',
                 bbox=myBBox,
                 size=(800, 800),
                 format='image/jpeg',
                 transparent=True
                 )
out = open('chateau.jpg', 'wb')
out.write(img.read())
out.close()

# + [markdown] id="btklwk1IfDUz" colab_type="text"
# # Traitements d'images avec PIL
#
# [Pillow](https://pillow.readthedocs.io/en/latest/handbook/tutorial.html) est une librairie qui permet de faire des traitements d'images. PIL est l'acronyme de Python Imaging Library.
#
#
#

# + [markdown] id="DJolzQgWujHW" colab_type="text"
# ## Charger une image en mémoire

# + id="gQ5EXxRhB9lQ" colab_type="code" colab={}
from PIL import Image
img = Image.open('chateau.jpg')

# + [markdown] id="EYouySz_uqx-" colab_type="text"
# ## Taille de l'image

# + id="QIbpZChiCM1M" colab_type="code" outputId="3f3413c7-315c-4c91-faba-99f444f72430" colab={"base_uri": "https://localhost:8080/", "height": 34}
img.size

# + [markdown] id="sokNHAP4unMn" colab_type="text"
# ## Visualiser l'image

# + id="Bb8b3XbAF-8U" colab_type="code" outputId="9090347b-ea28-4894-90de-327f20ce10e2" colab={"base_uri": "https://localhost:8080/", "height": 286}
# %matplotlib inline
from matplotlib.pyplot import imshow 
import numpy as np
imshow(img)

# + [markdown] id="5pr467jduwwR" colab_type="text"
# ## Transformer le RGB en HSV

# + [markdown] id="Ftbj_orvgeQL" colab_type="text"
# HSV=[Hue Saturation Value](https://fr.wikipedia.org/wiki/Teinte_Saturation_Valeur)
# - Hue:
#     - Couleur de l'arc-en-ciel
#     - Du rouge (Hue=0) vers le violet
#     - L'infrarouge et l'ultraviolet n'ont pas leur place dans cette échelle -> uniquement pour les couleurs visibles
#     - $H \in [0°,180°[$
# - Saturation:
#      - Intensité
#      - $S \in [0,1]$
#      - D'autant plus grand que l'on s'écarte du blanc, qui est la couleur qui intègre toutes les aures couleurs
#      - $\frac{max(R,G,B)-min(R,G,B)}{max(R,G,B)}$
#      - opencv --> 0 jusque 255
# - Valeur:
#     - Brillance (0=noir)
#     - $max(R,G,B)$
#     - $S \in [0,256]$
#     - opencv --> 0 jusque 255
#     
# Voir aussi articles:
#    - [Colour Rendering of Spectra](http://www.fourmilab.ch/documents/specrend/)
#    - [Hue, Value, Saturation](http://learn.leighcotnoir.com/artspeak/elements-color/hue-value-saturation/)

# + id="wMUdUyJ6oZ25" colab_type="code" outputId="3f7fbea1-8e49-4f09-d42b-7e5c7faddbb4" colab={"base_uri": "https://localhost:8080/", "height": 286}
import cv2
import numpy as np

## Read
img = cv2.imread("chateau.jpg")

## convert to hsv
hsv = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)

## Vérifier que la transformation inverse fonctionne correctement
imshow(cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB))

# + [markdown] id="FNqrmkmEu5cc" colab_type="text"
# ## Visualiser la distribution de HUE avec matplotlib

# + id="fUWd7pcNSs8v" colab_type="code" outputId="eceb4d7c-76bd-4ed5-f9b1-3b213463dd02" colab={"base_uri": "https://localhost:8080/", "height": 286}
from matplotlib.pyplot import imshow 

hue=hsv.copy()
hue[:,:,1]=255
hue[:,:,2]=255

imshow(cv2.cvtColor(hue, cv2.COLOR_HSV2RGB))

# + id="ElNXdarFHyah" colab_type="code" outputId="42bbae6a-87d3-49fb-c7ca-2c2239317f75" colab={"base_uri": "https://localhost:8080/", "height": 265}
import matplotlib.pyplot as plt

hue=hsv[:,:,0].flatten().tolist()

n, bins, patches = plt.hist(hue, 50)

bin_centers = (bins[:-1] + bins[1:])/2

col = bin_centers/180

cm = plt.cm.get_cmap('hsv')
#cm = plt.cm.get_cmap('gist_rainbow')
#cm = plt.cm.hsv

for c, p in zip(col, patches):
    plt.setp(p, 'facecolor', cm(c))
    
plt.show()

# + [markdown] id="UVBP4R1oRpGP" colab_type="text"
# ## Sélectionner portion d'image sur base de HUE
#
# Sélectionner le chemin en dolomie (bleu).

# + id="Lgs6xpAixCdA" colab_type="code" outputId="381095b1-4942-48a5-a609-dd7d16ac0f43" colab={"base_uri": "https://localhost:8080/", "height": 286}
seuil1=105
seuil2=130
image=hsv.copy()
image[:,:,0]=(((hsv[:,:,0]<seuil2)*1+(hsv[:,:,0]>seuil1)*1)==2)*75
image[:,:,1]=255
image[:,:,2]=255
imshow(cv2.cvtColor(image, cv2.COLOR_HSV2RGB))

# + [markdown] id="YzBzzLHdvGtf" colab_type="text"
# ## Matrice d'histogrammes pour HUE, SATURATION, VALUE

# + id="dzc3LxZ7pEYw" colab_type="code" outputId="ea607212-42f9-4852-cf15-a5cae75bcce6" colab={"base_uri": "https://localhost:8080/", "height": 729}
titles=['Hue','Saturation','Value']
f,a = plt.subplots(3,1,figsize=(10, 10))
a = a.ravel()
for idx,ax in enumerate(a):
    ax.hist(hsv[:,:,idx].flatten().tolist(),50)
    ax.set_title(titles[idx])
    ax.set_xlabel(titles[idx])
    ax.set_ylabel('Frequence')
plt.tight_layout()


# + [markdown] id="lxVI4DmT8NeD" colab_type="text"
#  # Transformer image en DataFrame

# + id="Mjvpgl9ghe7H" colab_type="code" outputId="e3dcdf25-127f-466a-a906-51e663ba1300" colab={"base_uri": "https://localhost:8080/", "height": 206}
import pandas as pd
myImageDF=pd.DataFrame({'Hue':       hsv[:,:,0].flatten().tolist(),
                        'Saturation':hsv[:,:,1].flatten().tolist(), 
                        'Value':     hsv[:,:,2].flatten().tolist()})
myImageDF.head()

# + [markdown] id="enYRjBJSgC3v" colab_type="text"
# # Clustering

# + [markdown] id="og5JPq9S8VxX" colab_type="text"
# ## Calculs

# + id="Ym0LuyhYfV6w" colab_type="code" colab={}
countOfClusters=4

# + id="9S02uG5Z7rhy" colab_type="code" colab={}
from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=countOfClusters)
kmeans.fit(myImageDF)

clusters = kmeans.predict(myImageDF)

# + [markdown] id="kx_-8C4X8bU7" colab_type="text"
# ## Ajouter une colonne au DataFrame avec le cluster de chaque point

# + id="nZP0d_cm71Uo" colab_type="code" outputId="3bcc67c7-1728-4745-9866-e602be15542c" colab={"base_uri": "https://localhost:8080/", "height": 206}
myImageDF['Cluster']=clusters
myImageDF.head()

# + id="nu-vYgNlh5JF" colab_type="code" outputId="f8b566dd-8c8c-4fc7-a7ae-940bbe8de455" colab={"base_uri": "https://localhost:8080/", "height": 286}
image=hsv.copy()
image[:,:,0]=clusters.reshape(800,800)/countOfClusters*180
image[:,:,1]=255
image[:,:,2]=255
imshow(cv2.cvtColor(image, cv2.COLOR_HSV2RGB))

# + [markdown] id="To0NCIUNRZUT" colab_type="text"
# ## Portion de territoire occupée par chaque couleur

# + id="s7efxRrx-tZJ" colab_type="code" outputId="ef22d4df-3159-4d75-f6c4-1ff94680253b" colab={"base_uri": "https://localhost:8080/", "height": 104}
myImageDF['Cluster'].value_counts()

# + id="d4_wzh3OPp6l" colab_type="code" outputId="35c6ed7f-8002-496e-af52-2b3be0e4573a" colab={"base_uri": "https://localhost:8080/", "height": 104}
myImageDF['Cluster'].value_counts()/myImageDF['Cluster'].count()

# + id="XCNVVtpWFqgX" colab_type="code" outputId="daf0b715-8945-4368-e2<PASSWORD>" colab={"base_uri": "https://localhost:8080/", "height": 265}
import matplotlib.pyplot as plt

n, bins, patches = plt.hist(clusters,bins=range(countOfClusters+1))

col = bins[:-1]/countOfClusters

cm = plt.cm.get_cmap('hsv')

for c, p in zip(col, patches):
    plt.setp(p, 'facecolor', cm(c))
    
plt.show()

StarcoderdataJupyterScriptsDedupFiltered