#!/usr/bin/env python
# coding: utf-8

# In[1]:


import fiona
import geopandas as gpd
import rasterio as rio
from rasterio.plot import show as rioshow
from rasterio import transform
from shapely.geometry import Point
from shapely.geometry import Polygon
from shapely.geometry import MultiPolygon
from shapely.geometry import LineString
from shapely.geometry import box
from shapely import wkt
from pyproj import Transformer

import matplotlib.pyplot as plt
import numpy as np
import os
import requests
from dotenv import load_dotenv

np.set_printoptions(suppress=True)


# In[2]:


# Enable KML support which is disabled by default
fiona.drvsupport.supported_drivers['libkml'] = 'rw'
fiona.drvsupport.supported_drivers['LIBKML'] = 'rw'


# ## Variables

# In[3]:


# Load environmaent variable
load_dotenv()


# In[4]:


geoCRS = 4326  # for output GeoJSON files (lon/lat)
projCRS = 3857  # for calculations (meters)
OPENTOPO_TOKEN = os.getenv('OPENTOPO_TOKEN')


# ## Initial Data Analysis

# ### Area of Interest

# In[5]:


# Store file path
kreminna_aoi_file = os.path.join('assets','geom','kreminna.kml')


# In[6]:


# Read and display the area of interest
kreminna_aoi = gpd.read_file(kreminna_aoi_file).to_crs(projCRS)
print(kreminna_aoi.crs)
kreminna_aoi


# In[7]:


kreminna_aoi.iloc[0,-1]


# In[8]:


# Remove z coordinate
kreminna_aoi['geometry'] = kreminna_aoi['geometry'].apply(lambda g: wkt.loads(wkt.dumps(g, output_dimension=2)))
kreminna_aoi


# In[9]:


# Remove unnecessary columns and rename one
kreminna_aoi = kreminna_aoi[['Name','geometry']].rename(columns={'Name':'name'})


# In[10]:


# Export area of interest
kreminna_aoi.to_crs(geoCRS).to_file('kreminna-clipped-aoi.geojson', driver='GeoJSON')


# ### Digital Elevation Model (DEM)

# In[11]:


# Box coordinates: west, south, east, north
# kreminna_aoi.bounds


# In[12]:


# Download Digital Elevation Model for the area (free API token at https://portal.opentopography.org/ is required)

# w, s, e, n = kreminna_aoi.to_crs(geoCRS).bounds.minx, kreminna_aoi.bounds.miny, kreminna_aoi.bounds.maxx, kreminna_aoi.bounds.maxy

# url = "https://portal.opentopography.org/API/globaldem"

# params = {
#     'demtype': 'SRTMGL1',
#     'south': s,
#     'north': n,
#     'west': w,
#     'east': e,
#     'outputFormat': 'GTiff',
#     'API_Key': OPENTOPO_TOKEN
# }

# response = requests.get(url, params=params)

# with open ('kreminna-DEM.tif', 'wb') as f:
#     f.write(response.content)


# In[13]:


# Store file path
kreminna_dem_file = 'kreminna-DEM.tif'


# In[14]:


# Read file contents
with rio.open(kreminna_dem_file) as f:
    print(f.meta)
    kreminna_dem_array = f.read(1)  # number of band (single here anyways)
    kreminna_dem_transform = f.transform
    kreminna_dem_crs = f.crs


# In[15]:


# Show dem and contours
rioshow(kreminna_dem_array, cmap='plasma', title='DEM'), rioshow(kreminna_dem_array, cmap='plasma', title='DEM Contour', contour=True)


# #### DEM to Contour in QGIS (+Smoothing)
# 
# The contour will be used later for the final visualization.
# 
# After importing the downloaded DEM the contours were modified. With the help of __"Generalizer"__ tool, most small objects present on the contour were removed (__"Remove small objects"__ algorithm). After that, the __"Boyle's Forward-Looking Algorithm"__ was used for line smoothing.  

# Before<br>
# <img src="assets/kreminna/kreminna-qgis-contour-raw.png" alt="w" width="700"/><br>
# 
# After<br>
# <img src="assets/kreminna/kreminna-qgis-contour-smooth.png" alt="w" width="700"/><br>

# In[16]:


# Store file path
kreminna_contour_file = os.path.join('assets','elevation','kreminna-elevation-contour.geojson')


# In[17]:


# Read and display data
kreminna_contour = gpd.read_file(kreminna_contour_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(kreminna_contour.crs)
kreminna_contour.plot()


# In[18]:


kreminna_contour.head()


# In[19]:


# Remove unnecessary columns and rename one
kreminna_contour = kreminna_contour[['ELEV','geometry']].rename(columns={'ELEV':'elevation'})


# In[20]:


# Check contours clip
kreminna_contour.clip(kreminna_aoi).plot()


# In[21]:


# Export clipped contour
kreminna_contour.clip(kreminna_aoi, keep_geom_type=True).to_crs(geoCRS).to_file('kreminna-clipped-contours.geojson', driver='GeoJSON')


# ### Occupied Area

# Some data cleaning and export as GeoJSON conducted using __[geojson.io](https://geojson.io/#map=2/0/20)__
# 
# __Before__: single file with different geometries
# 
# <img src="assets/kreminna/warzone-160624-geojson.io.png" alt="w" width="500"/><br>
# 
# __After__: 2 separate files, for areas (polygons) and fortifications (lines)
# 
# <img src="assets/kreminna/occupied-areas-160624-geojson.io.png" alt="oa" width="500"/><img src="assets/kreminna/fortification-lines-160624-geojson.io.png" alt="fl" width="500"/>

# In[22]:


# Store file path
warzone_output = os.path.join('assets','osint','occupied-areas-160624.geojson')


# In[23]:


# Read and display data
warzone_areas = gpd.read_file(warzone_output).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_areas.crs)
warzone_areas.plot()


# In[24]:


warzone_areas.shape


# In[25]:


warzone_areas


# In[26]:


# Remove z coordinate
warzone_areas['geometry'] = warzone_areas['geometry'].apply(lambda g: wkt.loads(wkt.dumps(g, output_dimension=2)))
warzone_areas.head(1)


# In[27]:


# Combine geometries
warzone_areas_d = gpd.geoseries.GeoSeries([geom for geom in warzone_areas.unary_union.geoms])
warzone_areas_d


# In[28]:


# Merge all polygons into a single multipolygon
warzone_area = warzone_areas.dissolve().replace('Luhansk Axis', 'occupied')
warzone_area


# In[29]:


# Display all occupied area
warzone_area.plot()


# In[30]:


# Display occupied area of interest
warzone_area.clip(kreminna_aoi).plot(color='r')


# In[31]:


# Export clipped occupied area
warzone_area.clip(kreminna_aoi, keep_geom_type=True).to_crs(geoCRS).to_file('kreminna-clipped-occupied.geojson', driver='GeoJSON')


# ### Trenches / Fortifications

# In[32]:


# Store file path
warzone_fortif_file= os.path.join('assets','osint','fortification-lines-160624.geojson')


# In[33]:


# Read and display data
warzone_fortif = gpd.read_file(warzone_fortif_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_fortif.crs)
warzone_fortif.plot()


# In[34]:


warzone_fortif.shape


# In[35]:


warzone_fortif


# In[36]:


# Check geometry types
warzone_fortif.geom_type.value_counts()


# In[37]:


# Remove unnecessary columns
warzone_fortif = warzone_fortif.drop(columns=['stroke','stroke-width','stroke-opacity','styleUrl'])


# In[38]:


# Check number of rows in 'unpacked' collections of geometries
warzone_fortif.explode(index_parts=False).shape


# In[39]:


# Store 'unpacked' geometries
warzone_fortif_decomposed = warzone_fortif.explode(ignore_index=True)
warzone_fortif_decomposed.plot()


# In[40]:


warzone_fortif_decomposed.tail()


# In[41]:


# Remove z coordianate
warzone_fortif_decomposed['geometry'] = warzone_fortif_decomposed['geometry'].apply(lambda g: wkt.loads(wkt.dumps(g, output_dimension=2)))
warzone_fortif_decomposed.head(1)


# In[42]:


# Check fortification types
warzone_fortif_decomposed['name'].value_counts()


# In[43]:


# Keep the fortification type name only
warzone_fortif_decomposed['name'] = warzone_fortif_decomposed['name'].apply(lambda x: x.split()[0])


# In[44]:


# Remove use of plural and singular forms of same word
warzone_fortif_decomposed['name'] = warzone_fortif_decomposed['name'].replace('Tankditch','Tankditches')


# In[45]:


# Plot all available fortifications data
fig, ax = plt.subplots(figsize=(6,8))
warzone_fortif_decomposed.plot(ax=ax, column='name', categorical=True, legend=True, legend_kwds={'loc':'center left', 'frameon':False})
plt.axis('off');


# In[46]:


# Plot fortifications data for the area of interest
warzone_fortif_decomposed.clip(kreminna_aoi).plot(column='name', categorical=True, legend=True, legend_kwds={'frameon':False, 'loc':(1,0.5)})
plt.axis('off');


# In[47]:


# Create a transformer to convert coordinates from EPSG:3857 to EPSG:4326
transformer = Transformer.from_crs(projCRS, geoCRS)


# In[48]:


# Convert the coordinates
x, y = transformer.transform(warzone_fortif_decomposed.centroid.x, warzone_fortif_decomposed.centroid.y)


# In[49]:


print('lat:', x[0])
print('lon:', y[0])


# In[50]:


# Extract lat and lon into new columns
warzone_fortif_decomposed['lat'] = x
warzone_fortif_decomposed['lon'] = y


# In[51]:


warzone_fortif_decomposed.head()


# In[52]:


# Check wsg84 coordinates
warzone_fortif_decomposed.to_crs(geoCRS).head(1)


# In[53]:


# Export clipped russian positions
warzone_fortif_decomposed.clip(kreminna_aoi, keep_geom_type=True).to_crs(geoCRS).to_file('kreminna-clipped-fortification-lines.geojson', driver='GeoJSON')


# ### Overpass Turbo Data

# Overpass Turbo uses the following format to define a bbox: __miny, minx, maxy, maxx__

# In[54]:


# Display area of interest bbox coordinates for overpass turbo query
kreminna_aoi.to_crs(geoCRS).bounds


# Therefore, bbox for current area of interest is: __(48.758341, 37.909804, 49.211731, 38.55444)__ 
# ***

# #### Roads

# __Overpass Turbo:__
# 
# ```
# [out:json][timeout:60];
# 
# way["highway"](48.758341,37.909804,49.211731,38.55444);
# 
# out geom;
# ```

# In[55]:


# Store file path
warzone_roads_file= os.path.join('assets','osint','kreminna-roads.geojson')


# In[56]:


# Read and display data
warzone_roads = gpd.read_file(warzone_roads_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_roads.crs)
warzone_roads.plot()


# In[57]:


warzone_roads.columns


# In[58]:


warzone_roads.iloc[0,:].values


# In[59]:


warzone_roads['highway'].value_counts(dropna=False)


# In[60]:


warzone_roads['highway'].isna().any()


# In[61]:


# Rename 'highway' column to a more generalized 'road'
warzone_roads = warzone_roads.rename(columns={'highway':'road'})


# In[62]:


# Replace None values in 'road' column with 'road' as their new value
warzone_roads['road'] = warzone_roads['road'].fillna('road')


# In[63]:


# Filter needed columns
warzone_roads = warzone_roads.loc[:,['id','road','geometry']]
warzone_roads.head(1)


# In[64]:


warzone_roads.isna().any()


# In[65]:


# Unpivot a GeoDataFrame
warzone_roads_unpivot = warzone_roads.melt(
    id_vars=['id','geometry'],
    value_vars='road', 
    value_name='sub_category', 
    var_name='category', 
    ignore_index=True
)

warzone_roads_unpivot.head()


# In[66]:


warzone_roads_unpivot.clip(kreminna_aoi).plot()


# In[67]:


# Export clipped roads
warzone_roads_unpivot.clip(kreminna_aoi).to_crs(geoCRS).to_file('kreminna-clipped-roads.geojson', driver='GeoJSON')


# #### Towns, villages, etc
# 
# __Overpass Turbo:__
# 
# ```
# [out:json][timeout:60];
# 
# node["place"]["place"!~"square"](48.758341,37.909804,49.211731,38.55444);
# 
# out geom;
# 
# ```

# In[68]:


# Store file path
warzone_places_file = os.path.join('assets','osint','kreminna-places.geojson')


# In[69]:


# Read and display data
warzone_places = gpd.read_file(warzone_places_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_places.crs)
warzone_places.plot()


# In[70]:


warzone_places.columns


# In[71]:


warzone_places.iloc[0,:].values


# In[72]:


warzone_places['name'].isna().value_counts()


# In[73]:


warzone_places['name:en'].isna().value_counts()


# In[74]:


# Filter needed columns
warzone_places_filtered = warzone_places.loc[:,['id','name:en','name','place','geometry']]
warzone_places_filtered.head(1)


# In[75]:


# Check which english place translations are missing
warzone_places_filtered[warzone_places_filtered['name:en'].isna()]


# In[76]:


# Fill missing place names
warzone_places_filtered.loc[38,'name:en'] = 'Vedmezhe'
warzone_places_filtered.loc[77,'name:en'] = 'Matroska'
warzone_places_filtered.loc[78,'name:en'] = 'Pivdenna'
warzone_places_filtered.loc[79,'name:en'] = 'Chykhyrove'
warzone_places_filtered[warzone_places_filtered['name:en'].isna()]


# In[77]:


warzone_places_filtered['place'].value_counts(dropna=False)


# In[78]:


# Replace None values in 'place' column with 'place' as their new value
warzone_places_filtered['place'] = warzone_places_filtered['place'].fillna('place')


# In[79]:


# Create a transformer to convert coordinates from EPSG:3857 to EPSG:4326
transformer = Transformer.from_crs(projCRS, geoCRS)

# Convert the coordinates
x, y = transformer.transform(warzone_places_filtered.centroid.x, warzone_places_filtered.centroid.y)


# In[80]:


# Extract lat and lon into new columns
warzone_places_filtered['lat'] = x
warzone_places_filtered['lon'] = y


# In[81]:


warzone_places_filtered.isna().any()


# In[82]:


# Unpivot GeoDataFrame
warzone_places_unpivot = warzone_places_filtered.melt(
    id_vars=['id','name:en','name','lon','lat','geometry'], 
    value_vars='place', 
    value_name='sub_category', 
    var_name='category', 
    ignore_index=True
)

warzone_places_unpivot.sample(10)


# In[83]:


warzone_places_unpivot.clip(kreminna_aoi).plot()


# In[84]:


# Export clipped places
warzone_places_unpivot.clip(kreminna_aoi).to_crs(geoCRS).to_file('kreminna-clipped-places.geojson', driver='GeoJSON')


# #### Buildings
# 
# ```
# [out:json][timeout:60];
# 
# nwr["building"](48.758341,37.909804,49.211731,38.55444);
# 
# out geom;
# ```

# In[85]:


# Store file path
warzone_buildings_file = os.path.join('assets','osint','kreminna-buildings.geojson')


# In[86]:


# Read and display data
warzone_buildings = gpd.read_file(warzone_buildings_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_buildings.crs)
warzone_buildings.plot()


# In[87]:


# Check geometry types
warzone_buildings.geom_type.value_counts(dropna=False)


# In[88]:


# Store point geometries instead of polygons
warzone_buildings['geometry'] = warzone_buildings.centroid
warzone_buildings.geom_type.value_counts(dropna=False)


# In[89]:


warzone_buildings.columns


# In[90]:


warzone_buildings.iloc[0,:].values


# In[91]:


warzone_buildings.isna().any()


# In[92]:


# Filter needed columns
warzone_buildings_filtered = warzone_buildings.loc[:,['id','geometry']]
warzone_buildings_filtered.head(1)


# In[93]:


# Create a transformer to convert coordinates from EPSG:3857 to EPSG:4326
transformer = Transformer.from_crs(projCRS, geoCRS)

# Convert the coordinates
x, y = transformer.transform(warzone_buildings_filtered.centroid.x, warzone_buildings_filtered.centroid.y)

# Extract lat and lon into new columns
warzone_buildings_filtered['lat'] = x
warzone_buildings_filtered['lon'] = y


# In[94]:


warzone_buildings_filtered.isna().any()


# In[95]:


warzone_buildings_filtered.clip(kreminna_aoi).plot()


# In[96]:


# Export clipped ponds
warzone_buildings_filtered.clip(kreminna_aoi).to_crs(geoCRS).to_file('kreminna-clipped-buildings.geojson', driver='GeoJSON')


# #### Lakes, ponds, etc
# 
# __Overpass Turbo:__
# 
# ```
# [out:json][timeout:60];
# 
# area["natural"="water"]["water"!~"river"](48.758341,37.909804,49.211731,38.55444);
# 
# out geom;
# 
# ```

# In[97]:


# Store file path
warzone_ponds_file = os.path.join('assets','osint','kreminna-ponds.geojson')


# In[98]:


# Read and display data
warzone_ponds = gpd.read_file(warzone_ponds_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_ponds.crs)
warzone_ponds.plot()


# In[99]:


warzone_ponds.columns


# In[100]:


warzone_ponds.iloc[0,:].values


# In[101]:


warzone_ponds.isna().any()


# In[102]:


warzone_ponds['water'].value_counts(dropna=False)


# In[103]:


# Replace None values in 'water' column with 'water' as their new value
warzone_ponds['water'] = warzone_ponds['water'].fillna('water')


# In[104]:


# Filter needed columns
warzone_ponds = warzone_ponds.loc[:,['id', '@id','water','geometry']]
warzone_ponds.head(1)


# In[105]:


warzone_ponds.isna().any()


# In[106]:


# Unpivot GeoDataFrame
warzone_ponds_unpivot = warzone_ponds.melt(
    id_vars=['id','geometry'],
    value_vars='water', 
    value_name='sub_category',
    var_name='category', 
    ignore_index=True
)

warzone_ponds_unpivot.sample(10)


# In[107]:


warzone_ponds_unpivot.clip(kreminna_aoi).plot()


# In[108]:


# Export clipped ponds
warzone_ponds_unpivot.clip(kreminna_aoi).to_crs(geoCRS).to_file('kreminna-clipped-ponds.geojson', driver='GeoJSON')


# #### Wetlands, swamps, etc
# 
# __Overpass Turbo:__
# 
# ```
# [out:json][timeout:60];
# 
# area["natural"="wetland"](48.758341,37.909804,49.211731,38.55444);
# 
# out geom;
# ```

# In[109]:


# Store file path
warzone_wetlands_file = os.path.join('assets','osint','kreminna-swamps.geojson')


# In[110]:


# Read and display data
warzone_wetlands = gpd.read_file(warzone_wetlands_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_wetlands.crs)
warzone_wetlands.plot()


# In[111]:


warzone_wetlands.columns


# In[112]:


warzone_wetlands.iloc[0,:].values


# In[113]:


warzone_wetlands.isna().any()


# In[114]:


warzone_wetlands.loc[:,'wetland'].value_counts(dropna=False)


# In[115]:


# Replace None values in 'wetland' column with 'wetland' as their new value
warzone_wetlands['wetland'] = warzone_wetlands['wetland'].fillna('wetland')


# In[116]:


# Filter needed columns
warzone_wetlands = warzone_wetlands.loc[:,['id', '@id','wetland','geometry']]
warzone_wetlands.head(1)


# In[117]:


warzone_wetlands.isna().any()


# In[118]:


# Unpivot GeoDataFrame
warzone_wetlands_unpivot = warzone_wetlands.melt(
    id_vars=['id', '@id','geometry'],
    value_vars='wetland', 
    value_name='sub_category', 
    var_name='category',
    ignore_index=True
)

warzone_wetlands_unpivot.head(10)


# In[119]:


warzone_wetlands_unpivot.clip(kreminna_aoi).plot()


# In[120]:


# Export clipped wetlands
warzone_wetlands_unpivot.clip(kreminna_aoi).to_crs(geoCRS).to_file('kreminna-clipped-wetlands.geojson', driver='GeoJSON')


# #### Rivers, streams
# 
# __Overpass Turbo:__
# 
# ```
# [out:json][timeout:60];
# 
# way["waterway"](48.758341,37.909804,49.211731,38.55444);
# 
# out geom;
# 
# ```

# In[121]:


# Store file path
warzone_rivers_file = os.path.join('assets','osint','kreminna-rivers.geojson')


# In[122]:


# Read and display data
warzone_rivers = gpd.read_file(warzone_rivers_file).to_crs(projCRS).dropna(how='all').reset_index(drop=True)
print(warzone_rivers.crs)
warzone_rivers.plot()


# In[123]:


warzone_rivers.columns


# In[124]:


warzone_rivers.iloc[0,:].values


# In[125]:


warzone_rivers.isna().any()


# In[126]:


warzone_rivers.loc[:,'waterway'].value_counts(dropna=False)


# In[127]:


warzone_rivers.loc[:,'width'].value_counts(dropna=False)


# In[128]:


# Replace None values in 'waterway' column with 'waterway' as their new value
warzone_rivers['waterway'] = warzone_rivers['waterway'].fillna('waterway')


# In[129]:


# Filter needed columns
warzone_rivers = warzone_rivers.loc[:,['id', '@id','waterway','geometry']]
warzone_rivers.head(1)


# In[130]:


# Unpivot GeoDataFrame
warzone_rivers_unpivot = warzone_rivers.melt(
    id_vars=['id', '@id','geometry'],
    value_vars='waterway', 
    value_name='sub_category', 
    var_name='category', 
    ignore_index=True
)

warzone_rivers_unpivot.sample(10)


# In[131]:


warzone_rivers_unpivot.clip(kreminna_aoi).plot()


# In[132]:


# Export clipped rivers
warzone_rivers_unpivot.clip(kreminna_aoi).to_crs(geoCRS).to_file('kreminna-clipped-rivers.geojson', driver='GeoJSON')