%matplotlib inline

import inspect, sys
import warnings; warnings.simplefilter('ignore')

# check pydov path
import pydov

"""
Explanation: Example of DOV search methods for soil data (bodemgegevens)

Use cases explained below

Introduction to the bodem-objects
Get bodemsites in a bounding box
Get bodemlocaties with specific properties
Get all direct and indirect bodemobservaties linked to a bodemlocatie
Get all bodemobservaties in a bodemmonster
Find all bodemlocaties where observations exist for organic carbon percentage in East-Flanders between 0 and 30 cm deep
Calculate carbon stock in Ghent in the layer 0 - 23 cm
End of explanation
"""


from pydov.search.bodemlocatie import BodemlocatieSearch
bodemlocatie = BodemlocatieSearch()

"""
Explanation: Get information about the datatype 'Bodemlocatie'
Other datatypes are also possible:
* Bodemsite: BodemsiteSearch
* Bodemmonster: BodemmonsterSearch
* Bodemobservatie: BodemobservatieSearch
End of explanation
"""


bodemlocatie.get_description()

"""
Explanation: A description is provided for the 'Bodemlocatie' datatype:
End of explanation
"""


fields = bodemlocatie.get_fields()

# print available fields
for f in fields.values():
    print(f['name'])

"""
Explanation: The different fields that are available for objects of the 'Bodemlocatie' datatype can be requested with the get_fields() method:
End of explanation
"""


fields['type']

"""
Explanation: You can get more information of a field by requesting it from the fields dictionary:
* name: name of the field
* definition: definition of this field
* cost: currently this is either 1 or 10, depending on the datasource of the field. It is an indication of the expected time it will take to retrieve this field in the output dataframe.
* notnull: whether the field is mandatory or not
* type: datatype of the values of this field
End of explanation
"""


fields['type']['values']

"""
Explanation: Optionally, if the values of the field have a specific domain the possible values are listed as values:
End of explanation
"""


from pydov.search.bodemsite import BodemsiteSearch
bodemsite = BodemsiteSearch()

from pydov.util.location import Within, Box

df = bodemsite.search(location=Within(Box(148000, 160800, 160000, 169500)))
df.head()

"""
Explanation: Example use cases
Get bodemsites in a bounding box
Get data for all the bodemsites that are geographically located completely within the bounds of the specified box.
The coordinates are in the Belgian Lambert72 (EPSG:31370) coordinate system and are given in the order of lower left x, lower left y, upper right x, upper right y.
The same methods can be used for other bodem objects.
End of explanation
"""


for pkey_bodemsite in set(df.pkey_bodemsite):
    print(pkey_bodemsite)

"""
Explanation: The dataframe contains a list of bodemsites. The available data are flattened to represent unique attributes per row of the dataframe.
Using the pkey_bodemsite field one can request the details of this bodemsite in a webbrowser:
End of explanation
"""


[i for i,j in inspect.getmembers(sys.modules['owslib.fes'], inspect.isclass) if 'Property' in i]

"""
Explanation: Get bodemlocaties with specific properties
Next to querying bodem objects based on their geographic location within a bounding box, we can also search for bodem objects matching a specific set of properties. 
The same methods can be used for all bodem objects.
For this we can build a query using a combination of the 'Bodemlocatie' fields and operators provided by the WFS protocol.
A list of possible operators can be found below:
End of explanation
"""


from owslib.fes import PropertyIsEqualTo

query = PropertyIsEqualTo(propertyname='bodemstreek',
                          literal='Zandstreek')
df = bodemlocatie.search(query=query, max_features=10)

df.head()

"""
Explanation: In this example we build a query using the PropertyIsEqualTo operator to find all bodemlocaties with bodemstreek 'zandstreek'.
We use max_features=10 to limit the results to 10.
End of explanation
"""


for pkey_bodemlocatie in set(df.pkey_bodemlocatie):
    print(pkey_bodemlocatie)

"""
Explanation: Once again we can use the pkey_bodemlocatie as a permanent link to the information of these bodemlocaties:
End of explanation
"""


from pydov.search.bodemobservatie import BodemobservatieSearch
from pydov.search.bodemlocatie import BodemlocatieSearch
bodemobservatie = BodemobservatieSearch()
bodemlocatie = BodemlocatieSearch()

from owslib.fes import PropertyIsEqualTo
from pydov.util.query import Join

bodemlocaties = bodemlocatie.search(query=PropertyIsEqualTo(propertyname='naam', literal='VMM_INF_52'), 
                                   return_fields=('pkey_bodemlocatie',))

bodemobservaties = bodemobservatie.search(query=Join(bodemlocaties, 'pkey_bodemlocatie'))
bodemobservaties.head()

"""
Explanation: Get all direct and indirect bodemobservaties in bodemlocatie
Get all bodemobservaties in a specific bodemlocatie.
Direct means bodemobservaties directly linked with a bodemlocatie.
Indirect means bodemobservaties linked with child-objects of the bodemlocatie, like bodemmonsters.
End of explanation
"""


from pydov.search.bodemmonster import BodemmonsterSearch
bodemmonster = BodemmonsterSearch()

bodemmonsters = bodemmonster.search(query=PropertyIsEqualTo(propertyname = 'identificatie', literal='A0057359'),
                                   return_fields=('pkey_bodemmonster',))

bodemobservaties = bodemobservatie.search(query=Join(bodemmonsters, on = 'pkey_parent', using='pkey_bodemmonster'))
bodemobservaties.head()

"""
Explanation: Get all bodemobservaties in a bodemmonster
Get all bodemobservaties linked with a bodemmonster
End of explanation
"""


from owslib.fes import PropertyIsEqualTo
from pydov.util.query import Join

from pydov.search.bodemclassificatie import BodemclassificatieSearch
from pydov.search.bodemlocatie import BodemlocatieSearch

bodemclassificatie = BodemclassificatieSearch()
bl_Scbz = bodemclassificatie.search(query=PropertyIsEqualTo('bodemtype', 'Scbz'), return_fields=['pkey_bodemlocatie'])

bodemlocatie = BodemlocatieSearch()
bl = bodemlocatie.search(query=Join(bl_Scbz, 'pkey_bodemlocatie'))
bl.head()

"""
Explanation: Find all soil locations with a given soil classification
Get all soil locations with a given soil classification:
End of explanation
"""


from pydov.search.bodemobservatie import BodemobservatieSearch

bodemobservatie = BodemobservatieSearch()
obs = bodemobservatie.search(query=Join(bl_Scbz, 'pkey_bodemlocatie'), max_features=10)
obs.head()

"""
Explanation: We can also get their observations:
End of explanation
"""


from pydov.search.bodemlocatie import BodemlocatieSearch
from pydov.search.bodemdiepteinterval import BodemdiepteintervalSearch
from pydov.util.query import Join
from owslib.fes import PropertyIsEqualTo

bodemlocatie = BodemlocatieSearch()
bodemdiepteinterval = BodemdiepteintervalSearch()

bodemlocaties = bodemlocatie.search(query=PropertyIsEqualTo(propertyname='naam', literal='VMM_INF_52'),
                                    return_fields=('pkey_bodemlocatie',))

bodemdiepteintervallen = bodemdiepteinterval.search(
    query=Join(bodemlocaties, on='pkey_bodemlocatie'))
bodemdiepteintervallen

"""
Explanation: Get all depth intervals and observations from a soil location
End of explanation
"""


from pydov.search.bodemobservatie import BodemobservatieSearch

bodemobservatie = BodemobservatieSearch()

bodemobservaties = bodemobservatie.search(query=Join(
    bodemdiepteintervallen, on='pkey_parent', using='pkey_diepteinterval'))
bodemobservaties.head()

"""
Explanation: And get their observations:
End of explanation
"""


from owslib.etree import etree
from owslib.wfs import WebFeatureService
from pydov.util.location import (
    GmlFilter,
    Within,
)

provinciegrenzen = WebFeatureService(
    'https://geoservices.informatievlaanderen.be/overdrachtdiensten/VRBG/wfs',
    version='1.1.0')

provincie_filter = PropertyIsEqualTo(propertyname='NAAM', literal='Oost-Vlaanderen')
provincie_poly = provinciegrenzen.getfeature(
    typename='VRBG:Refprv',
    filter=etree.tostring(provincie_filter.toXML()).decode("utf8")).read()

"""
Explanation: Find all bodemlocaties where observations exist for organic carbon percentage in East-Flanders between 0 and 30 cm deep
Get boundaries of East-Flanders by using a WFS
End of explanation
"""


from owslib.fes import PropertyIsEqualTo
from owslib.fes import And

from pydov.search.bodemobservatie import BodemobservatieSearch

bodemobservatie = BodemobservatieSearch()

# Select only layers with the boundaries 10-30
bodemobservaties = bodemobservatie.search(
        location=GmlFilter(provincie_poly, Within),
        query=And([
            PropertyIsEqualTo(propertyname="parameter", literal="Organische C - percentage"),
            PropertyIsEqualTo(propertyname="diepte_tot_cm", literal = '30'),
            PropertyIsEqualTo(propertyname="diepte_van_cm", literal = '0')
        ]))


bodemobservaties.head()

"""
Explanation: Get bodemobservaties in East-Flanders with the requested properties
End of explanation
"""


from pydov.search.bodemlocatie import BodemlocatieSearch
from pydov.util.query import Join
import pandas as pd

# Find bodemlocatie information for all observations
bodemlocatie = BodemlocatieSearch()
bodemlocaties = bodemlocatie.search(query=Join(bodemobservaties, on = 'pkey_bodemlocatie', using='pkey_bodemlocatie'))

# remove x, y, mv_mtaw from observatie dataframe to prevent duplicates while merging
bodemobservaties = bodemobservaties.drop(['x', 'y', 'mv_mtaw'], axis=1)

# Merge the bodemlocatie information together with the observation information
merged = pd.merge(bodemobservaties, bodemlocaties, on="pkey_bodemlocatie", how='left')

merged.head()

"""
Explanation: Now we have all observations with the requested properties. 
Next we need to link them with the bodemlocatie
End of explanation
"""


import folium
from folium.plugins import MarkerCluster
from pyproj import Transformer

# convert the coordinates to lat/lon for folium
def convert_latlon(x1, y1):
    transformer = Transformer.from_crs("epsg:31370", "epsg:4326", always_xy=True)
    x2,y2 = transformer.transform(x1, y1)
    return x2, y2

#convert coordinates to wgs84
merged['lon'], merged['lat'] = zip(*map(convert_latlon, merged['x'], merged['y']))

# Get only location and value
loclist = merged[['lat', 'lon']].values.tolist()

# initialize the Folium map on the centre of the selected locations, play with the zoom until ok
fmap = folium.Map(location=[merged['lat'].mean(), merged['lon'].mean()], zoom_start=10)
marker_cluster = MarkerCluster().add_to(fmap)
for loc in range(0, len(loclist)):
    popup = 'Bodemlocatie: ' + merged['pkey_bodemlocatie'][loc] 
    popup = popup + '<br> Bodemobservatie: ' + merged['pkey_bodemobservatie'][loc]
    popup = popup + '<br> Value: ' + merged['waarde'][loc] + "%"
    folium.Marker(loclist[loc], popup=popup).add_to(marker_cluster)
fmap

"""
Explanation: To export the results to CSV, you can use for example: 
python
merged.to_csv("test.csv")
We can plot also the results on a map
This can take some time!
End of explanation
"""


from owslib.etree import etree
from owslib.fes import PropertyIsEqualTo
from owslib.wfs import WebFeatureService
from pydov.util.location import (
    GmlFilter,
    Within,
)

stadsgrenzen = WebFeatureService(
    'https://geoservices.informatievlaanderen.be/overdrachtdiensten/VRBG/wfs',
    version='1.1.0')

gent_filter = PropertyIsEqualTo(propertyname='NAAM', literal='Gent')
gent_poly = stadsgrenzen.getfeature(
    typename='VRBG:Refgem',
    filter=etree.tostring(gent_filter.toXML()).decode("utf8")).read()


"""
Explanation: Calculate carbon stock in Ghent in the layer 0 - 23 cm
At the moment, there are no bulkdensities available. As soon as there are observations with bulkdensities, this example can be used to calculate a carbon stock in a layer.
Get boundaries of Ghent using WFS
End of explanation
"""


from owslib.fes import PropertyIsEqualTo, PropertyIsGreaterThan, PropertyIsLessThan
from owslib.fes import And

from pydov.search.bodemobservatie import BodemobservatieSearch

bodemobservatie = BodemobservatieSearch()

# all layers intersect the layer 0-23cm
carbon_observaties = bodemobservatie.search(
        location=GmlFilter(gent_poly, Within),
        query=And([
            PropertyIsEqualTo(propertyname="parameter", literal="Organische C - percentage"),
            PropertyIsGreaterThan(propertyname="diepte_tot_cm", literal = '0'),
            PropertyIsLessThan(propertyname="diepte_van_cm", literal = '23')
        ]),
        return_fields=('pkey_bodemlocatie', 'waarde'))
carbon_observaties = carbon_observaties.rename(columns={"waarde": "organic_c_percentage"})
carbon_observaties.head()


"""
Explanation: First get all observations in Ghent for organisch C percentage in requested layer
End of explanation
"""


density_observaties = bodemobservatie.search(
        location=GmlFilter(gent_poly, Within),
        query=And([
            PropertyIsEqualTo(propertyname="parameter", literal="Bulkdensiteit - gemeten"),
            PropertyIsGreaterThan(propertyname="diepte_tot_cm", literal = '0'),
            PropertyIsLessThan(propertyname="diepte_van_cm", literal = '23')
        ]),
        return_fields=('pkey_bodemlocatie', 'waarde'))

density_observaties = density_observaties.rename(columns={"waarde": "bulkdensity"})
density_observaties.head()

"""
Explanation: Then get all observations in Ghent for bulkdensity in requested layer
End of explanation
"""


import pandas as pd

merged = pd.merge(carbon_observaties, density_observaties, on="pkey_bodemlocatie")

merged.head()

"""
Explanation: Merge results together based on their bodemlocatie. Only remains the records where both parameters exists
End of explanation
"""


from owslib.fes import PropertyIsLike

query = PropertyIsLike(propertyname='naam',
                       literal='KART_PROF_%', wildCard='%')
df = bodemlocatie.search(query=query, max_features=10)

df.head()

"""
Explanation: Filter Aardewerk soil locations
Since we know that Aardewerk soil locations make use of a specific suffix, a query could be built filtering these out.
Since we only need to match a partial string in the name, we will build a query using the PropertyIsLike operator to find all Aardewerk bodemlocaties.
We use max_features=10 to limit the results to 10.
End of explanation
"""


for pkey_bodemlocatie in set(df.pkey_bodemlocatie):
    print(pkey_bodemlocatie)

"""
Explanation: As seen in the soil data example, we can use the pkey_bodemlocatie as a permanent link to the information of these bodemlocaties:
End of explanation
"""


