Temporal variability of climate time Series using MODIS satellite data
Introduction
The colab is avaiable here
Install planetary_computer and import other libraries
# Install libraries
!pip install planetary_computer
!pip install odc-stac
!pip install osmnx
!pip install rioxarray
!pip install zarr
!pip install geojson
# Import libraries
from urllib.parse import urlparse
from pystac.extensions.eo import EOExtension as eo
import matplotlib.pyplot as plt
import planetary_computer
import pystac
import pystac_client
import odc.stac
import rich.table
import osmnx as ox
import folium
import numpy as np
import rioxarray
import rasterio
import xarray as xr
import zarr
import os
Use pystac to import a stac catalog from microsoft planetary computer
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
Choose an area, time of interest and dataset from planetary computer
bbox = [3.771744,43.610725,4.335823,43.814463] # 738785.0528,6274273.2065,779106.8653,6298491.2641
year = "2021"
months = {
"January": "01",
"February": "02",
"March": "03",
"April": "04",
"May": "05",
"June": "06",
"July": "07",
"August": "08",
"September": "09",
"October": "10",
"november": "11",
"December": "12",
}
items = dict()
# Fetch the collection of interest and print available items
for name, number in months.items():
datetime = f"{year}-{number}"
search = catalog.search(
collections=["modis-11A1-061"],
bbox=bbox,
datetime=datetime,
)
items[name] = search.item_collection()[0]
print(items)
print(len(items))
{'January': <Item id=MYD11A1.A2021031.h18v04.061.2021041154946>, 'February': <Item id=MYD11A1.A2021059.h18v04.061.2021062101022>, 'March': <Item id=MYD11A1.A2021090.h18v04.061.2021092202346>, 'April': <Item id=MYD11A1.A2021120.h18v04.061.2021122015755>, 'May': <Item id=MYD11A1.A2021151.h18v04.061.2021152213240>, 'June': <Item id=MYD11A1.A2021181.h18v04.061.2021183192128>, 'July': <Item id=MYD11A1.A2021212.h18v04.061.2021213191344>, 'August': <Item id=MYD11A1.A2021243.h18v04.061.2021244183922>, 'September': <Item id=MYD11A1.A2021273.h18v04.061.2021277221831>, 'October': <Item id=MYD11A1.A2021304.h18v04.061.2021306213324>, 'november': <Item id=MYD11A1.A2021334.h18v04.061.2021335171453>, 'December': <Item id=MYD11A1.A2021365.h18v04.061.2022004015033>}
12
Choose the right asset
Each item has several available assets, including the original HDF file and a Cloud-optimized GeoTIFF of each subdataset. In our case, our interest is the ` LST_Day_1km` asset, that corresponds to the temperature mesured during daytime
rich_table = rich.table.Table("Key", "Title")
for key, asset in items["March"].assets.items():
rich_table.add_row(key, asset.title)
rich_table
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Key ┃ Title ┃ ┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ hdf │ Source data containing all bands │ │ QC_Day │ Quality control for daytime LST and emissivity │ │ Emis_31 │ Band 31 emissivity │ │ Emis_32 │ Band 32 emissivity │ │ QC_Night │ Quality control for nighttime LST and emissivity │ │ metadata │ Federal Geographic Data Committee (FGDC) Metadata │ │ LST_Day_1km │ Daily daytime 1km grid Land-surface Temperature │ │ Clear_day_cov │ Day clear-sky coverage │ │ Day_view_angl │ View zenith angle of daytime Landsurface Temperature │ │ Day_view_time │ (local solar) Time of daytime Land-surface Temperature observation │ │ LST_Night_1km │ Daily nighttime 1km grid Land-surface Temperature │ │ Clear_night_cov │ Night clear-sky coverage │ │ Night_view_angl │ View zenith angle of nighttime Landsurface Temperature │ │ Night_view_time │ (local solar) Time of nighttime Landsurface Temperature observation │ │ tilejson │ TileJSON with default rendering │ │ rendered_preview │ Rendered preview │ └──────────────────┴─────────────────────────────────────────────────────────────────────┘
The MODIS coordinate reference system is a sinusoidal grid, which means that views in a naïve XY raster look skewed. For visualization purposes, we reproject to a spherical Mercator projection for intuitive, north-up visualization.
data = odc.stac.load(
items.values(),
crs="EPSG:2154",
bands="LST_Day_1km",
resolution=500,
bbox=bbox,
)
# Get data data of July as a sample
raster = items["July"].assets["LST_Day_1km"].extra_fields["raster:bands"]
# Use raster scale to redefine temperature in celsius
data = (data["LST_Day_1km"] * raster[0]["scale"]) - 273.15
data
<xarray.DataArray 'LST_Day_1km' (time: 12, y: 48, x: 92)>
array([[[-273.15, -273.15, -273.15, ..., 16.55, 16.71, 16.71],
[-273.15, -273.15, -273.15, ..., 16.55, 16.71, 16.71],
[-273.15, -273.15, -273.15, ..., 16.73, 16.69, 16.69],
...,
[ 15.51, 15.21, 15.21, ..., 11.71, 11.35, 11.35],
[ 15.47, 15.55, 15.55, ..., 11.71, 11.35, 11.35],
[ 15.47, 15.55, 15.55, ..., 13.47, 14.21, 14.21]],
[[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15],
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15],
[-273.15, -273.15, -273.15, ..., -273.15, 20.53, 20.53],
...,
[ 16.67, 17.83, 17.83, ..., 13.59, 14.23, 14.23],
[-273.15, 17.83, 17.83, ..., 13.59, 14.23, 14.23],
[-273.15, 17.83, 17.83, ..., 15.63, 15.97, 15.97]],
[[ 23.83, 23.27, 23.27, ..., 31.33, 31.65, 31.65],
[ 24.17, 23.93, 23.93, ..., 31.33, 31.65, 31.65],
[ 24.17, 23.93, 23.93, ..., 31.13, 30.55, 30.55],
...,
...
...,
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15],
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15],
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15]],
[[ 8.65, 8.65, 8.65, ..., 9.33, 9.35, 9.35],
[ 8.63, 8.63, 8.63, ..., 9.33, 9.35, 9.35],
[ 8.63, 8.63, 8.63, ..., 9.55, 9.47, 9.47],
...,
[ 10.17, 10.23, 10.23, ..., 6.01, 6.09, 6.09],
[ 10.21, 10.17, 10.17, ..., 6.01, 6.09, 6.09],
[ 10.21, 10.17, 10.17, ..., 6.91, 7.67, 7.67]],
[[ 16.45, 16.07, 16.07, ..., -273.15, -273.15, -273.15],
[ 16.59, 16.17, 16.17, ..., -273.15, -273.15, -273.15],
[ 16.59, 16.17, 16.17, ..., -273.15, -273.15, -273.15],
...,
[ 16.65, 16.99, 16.99, ..., -273.15, -273.15, -273.15],
[ 16.63, 16.79, 16.79, ..., -273.15, -273.15, -273.15],
[ 16.63, 16.79, 16.79, ..., -273.15, -273.15, -273.15]]])
Coordinates:
* y (y) float64 6.303e+06 6.302e+06 ... 6.28e+06 6.279e+06
* x (x) float64 7.622e+05 7.628e+05 ... 8.072e+05 8.078e+05
spatial_ref int32 2154
* time (time) datetime64[ns] 2021-01-31 2021-02-28 ... 2021-12-31Interpolate data and fill missing values
There are many nodata values. Lets interpolate missing values using rioxarray. But me must first define our nodata value. In this case, we can define as -273.15, since Modis defines nodata value as 0
data_with_nodata = data.rio.write_nodata(-273.15, inplace=True)
print(f'No data value is {data_with_nodata.rio.nodata}')
No data value is -273.15
We plot the month of mars
data_with_nodata[2].plot()
<matplotlib.collections.QuadMesh at 0x798254418730>
With the nodata defined, we can interpolate na’s using the nearest neighbour method and fill missing data
# Interpolate missing values
data_interpolated=data_with_nodata.rio.interpolate_na(method='nearest')
# Plot march again
data_interpolated[2].plot()
<matplotlib.collections.QuadMesh at 0x798254a2e380>
Plot all the months together
g = data_interpolated.plot.imshow(cmap="Spectral_r", col="time", vmin=-10, vmax=60, size=4,col_wrap=4)
datetimes = data.time.to_pandas().dt.strftime("%B")
for ax, datetime in zip(g.axes.flat, datetimes):
ax.set_title(datetime)
<ipython-input-39-d1854156f2f9>:4: DeprecationWarning: self.axes is deprecated since 2022.11 in order to align with matplotlibs plt.subplots, use self.axs instead.
for ax, datetime in zip(g.axes.flat, datetimes):
Plot images with folium
Create bbox with Fiona
import fiona
import rasterio
import rasterio.mask
import json
import geojson
from fiona.crs import from_epsg
from shapely.geometry import mapping
import json
import requests
response = requests.get('http://app3.ouicodedata.com/collections/agropolis-zone/items')
with open('./response.json', 'w') as outfile:
json.dump(response.json(), outfile)
# Load JSON data from a file
with open('./response.json','r') as infile:
with open('./agropoilis_zone.geojson', 'w') as outfile:
data = json.load(infile)
geojson.dump(data,outfile)
with fiona.open("./agropoilis_zone.geojson", "r",crs=from_epsg(2154)) as geojson:
shapes = [feature["geometry"] for feature in geojson]
print(from_epsg(2154))
EPSG:2154
import os
from rasterio import windows
import random
from rasterio import features
from rasterio import warp
from rasterio.windows import get_data_window
import rasterio.mask
for i,d in enumerate(data_interpolated,start=1):
if not os.path.exists('temp'): os.makedirs('temp')
temp_file_name = '_'.join(['./temp/modis_tile',list(months.values())[i-1],'.tiff'])
d.rio.to_raster(temp_file_name)
with rasterio.open(temp_file_name) as src:
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
out_meta = src.meta
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
if not os.path.exists('cropped'): os.makedirs('cropped')
croped_file_name = '_'.join(['./cropped/cropped_modis_tile',list(months.values())[i-1],'.tiff'])
with rasterio.open(croped_file_name, "w", **out_meta) as dest:
dest.write(out_image)
Export images and plot then together with Folium
from pyproj import Transformer
f = folium.Figure(width=700, height=500)
m = folium.Map(location=(43.7237837,4.0003508)).add_to(f)
from rasterio.windows import Window
for i,d in enumerate(data,start=1):
temp_file_name = '_'.join(['./temp/modis_tile',list(months.values())[i-1],'.tiff'])
dst_crs = 'EPSG:4326'
with rasterio.open(temp_file_name) as src:
img = src.read()
src_crs = src.crs['init'].upper()
min_lon, min_lat, max_lon, max_lat = src.bounds
## Conversion from UTM to WGS84 CRS
bounds_orig = [[min_lat, min_lon], [max_lat, max_lon]]
bounds_fin = []
for item in bounds_orig:
#converting to lat/lon
lat = item[0]
lon = item[1]
proj = Transformer.from_crs(int(src_crs.split(":")[1]), int(dst_crs.split(":")[1]), always_xy=True)
lon_n, lat_n = proj.transform(lon, lat)
bounds_fin.append([lat_n, lon_n])
# Finding the centre latitude & longitude
centre_lon = bounds_fin[0][1] + (bounds_fin[1][1] - bounds_fin[0][1])/2
centre_lat = bounds_fin[0][0] + (bounds_fin[1][0] - bounds_fin[0][0])/2
layer = folium.raster_layers.ImageOverlay(
image = img.transpose(1, 2, 0),
name=list(months.keys())[i-1], # add a name to the layer
bounds = bounds_fin,
opacity=0.7,
interactive=False,
cross_origin=False,
zindex=1,
)
# Overlay the image
m.add_child(layer)
folium.LayerControl().add_to(m)
m
<iframe srcdoc=”<!DOCTYPE html> <html> <head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8" />
<script>
L_NO_TOUCH = false;
L_DISABLE_3D = false;
</script>
<style>html, body {width: 100%;height: 100%;margin: 0;padding: 0;}</style>
<style>#map {position:absolute;top:0;bottom:0;right:0;left:0;}</style>
<script src="https://cdn.jsdelivr.net/npm/leaflet@1.9.3/dist/leaflet.js"></script>
<script src="https://code.jquery.com/jquery-1.12.4.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/bootstrap@5.2.2/dist/js/bootstrap.bundle.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/Leaflet.awesome-markers/2.0.2/leaflet.awesome-markers.js"></script>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/leaflet@1.9.3/dist/leaflet.css"/>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/bootstrap@5.2.2/dist/css/bootstrap.min.css"/>
<link rel="stylesheet" href="https://netdna.bootstrapcdn.com/bootstrap/3.0.0/css/bootstrap.min.css"/>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/@fortawesome/fontawesome-free@6.2.0/css/all.min.css"/>
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/Leaflet.awesome-markers/2.0.2/leaflet.awesome-markers.css"/>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/python-visualization/folium/folium/templates/leaflet.awesome.rotate.min.css"/>
<meta name="viewport" content="width=device-width,
initial-scale=1.0, maximum-scale=1.0, user-scalable=no" />
<style>
#map_c411caddeadbebf70f7cebfd52330cd5 {
position: relative;
width: 100.0%;
height: 100.0%;
left: 0.0%;
top: 0.0%;
}
.leaflet-container { font-size: 1rem; }
</style>
<style>
.leaflet-image-layer {
/* old android/safari*/
image-rendering: -webkit-optimize-contrast;
image-rendering: crisp-edges; /* safari */
image-rendering: pixelated; /* chrome */
image-rendering: -moz-crisp-edges; /* firefox */
image-rendering: -o-crisp-edges; /* opera */
-ms-interpolation-mode: nearest-neighbor; /* ie */
}
</style>
</head> <body>
<div class="folium-map" id="map_c411caddeadbebf70f7cebfd52330cd5" ></div>
</body> <script>
var map_c411caddeadbebf70f7cebfd52330cd5 = L.map(
"map_c411caddeadbebf70f7cebfd52330cd5",
{
center: [43.7237837, 4.0003508],
crs: L.CRS.EPSG3857,
zoom: 10,
zoomControl: true,
preferCanvas: false,
}
);
var tile_layer_8ebb7bdf73674872bc467e44a609f99a = L.tileLayer(
"https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",
{"attribution": "Data by \u0026copy; \u003ca target=\"_blank\" href=\"http://openstreetmap.org\"\u003eOpenStreetMap\u003c/a\u003e, under \u003ca target=\"_blank\" href=\"http://www.openstreetmap.org/copyright\"\u003eODbL\u003c/a\u003e.", "detectRetina": false, "maxNativeZoom": 18, "maxZoom": 18, "minZoom": 0, "noWrap": false, "opacity": 1, "subdomains": "abc", "tms": false}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_06ec3d8b1ee2472e0dee5d2551cc2776 = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAJmUlEQVR42uXbBY8USxSG4f7DBIIEJ0hwd0Jwd3d3d3d3CAGChzA330neztnaap3ZAXJv8iULd5nuerq66lR1TzJ58uRGUYYNG9aYMWOGZd26dZYtW7bkZsmSJZkZM2ZMmgkTJjSdlStX5ubIkSOW2P/bu3dvh6xfv96ydu3aNCNGjCidBQsWWG7cuGH5/Pmz5du3b5YkD3r8+PGGXQW8LHQrwfPQwW4GvCx6z549LW0BrwrdDvSq2HXB+/btawH8x48faYDOBR84cGCawYMHWyZNmmTZvHmz5dy5c5bz589b6kB3FfiaNWssYG/durUT9vLlyyuB56FXAp87d24mdjPgZbFbCQ7awYMHLceOHbOE4MKuCs4dHgPv3r17h3z48CHN+/fvO8TASYgdguvEycmTJy379++3LF682FIF+k+Ag10X3KOH0JXBY+iAe+ws8DrYrQLnPHwuX75sOXHiRFPgHlsZN26cZeTIkYXgjx8/tjx9+tTSCTxEj2Eru3btsujkxo4da/lT4DFswIVdFTyvd4fg5PDhw5YLFy5Yrl27lmIXgvfp0yfNokWLugScO4e0Evzo0aMWzS2UZ/v27bNs2rQpzenTp9Mw/OiiTJs2rUPqgvtcvXrVUgpcCSsZykPh8ztloHVyIXZZ8BUrVkQTw24nONiVwbt169YB24OH1UwIzkSlxJBJDPvWrVuWssB5AXvHjh2W69evW4BljFZOnTqVhqFn+/btjSlTpnSISkwPjkEMnOoNYP5MqZpQwgm7CDxEP3DggMVjE4+ch+3BlTrIPkAXgW/YsKETuLA9OEgkCxwHhc/juJTNVHG54IMGDWosXLjQwkGXLl1qUeOYaHbv3p1m4sSJtbG7Ap3Fz5UrVyycs8AVkHXuHjzE9ujCloEA6ZClwXXQ+fPndwIXdlVwYcfAq2A3C65xXMDNgCvTp09P44c52h+Cg86ccenSJQvH37lzpyVhLKZsoqgfMmSIZdasWRYmHW7JZcuWpatKTj4GXhW7GXQmTjWQnw8dOmRhTKXH8ffbtm2z8BnNgCtdAi7sMuB1sZsBp4Ex8Nu3b6cQtGfPnj2dPmf06NGZARvgefPmWfjzxo0bLVw47jTdVVooJshT5Gt4UdiIETYzraJbj4vEvgqfsXr1asvQoUMtIfSZM2csrQanpBMgDaVaYUwVdhlwAP8YuCLkZsGrYJcFB5qoUYAI2/+/sN5mElNi4Fnw+uzZs2d3ypw5cyzhEOMvrFa+SbilyRjN3zMhsndCT1KAphRk9clnhdg3b94sDV6EHmKH4B63CFwJoZVRo0ZZevXqlSaG3WXg/tYtAw402K0A9w8VfFh4sIABWA1VuNWZo3iwImw9QCkLnoXOqhR4TFQa4qbJOhk+fHhDYXLkxBguwNVwwskDytMeej+1KOOnoM+ePWupAu7LsyLoZsH9U6s64OE2wF8Dfu/ePUu4x+DjocnUqVMLsf2kzvEAp+FULYzFanP4mBBoPVpU+vfvbwnBPXoWuN9bYRGkDpeCs1xlWUt9KuDjx49baCC9jovDJMoFoZfpQj18+NCSBx6D9uCxxLDbCc7YHwMPN7OErba/efOmOXB6eliX6yJQjt2/f99y9+5di9+/yINWVMYVwSt+D4eGAs85c1tr1czEyq4hxwCNYYdqSwmhfTw4x6cMVehsBg40twI9AHC/rPcrOaUZ8CJsLcM54SroXQVOG2PgoMewX758mebTp0+dwf1+MQdj1ckQQjXCrchKE3C/18xkyVUugs4Dz0NnHQEwS3g6Dm3yFRlVFu1gXO7du7eF9utuLgJXwueXpAM4m1OrVq2ycGLa5/5T4Gw0xcCz4FsFzsJP4L4c5bFcFrbG6SxwXgb6/v17IwEv3NsWOLfYgAEDLOwHZy1jNd7TYMDBo+Tjdg43+YnmihD8wYMHheDcobSDyomhRaDgcA7h2gNsUhU8C/2vBWdyjoEXoeeBq1LQ4g08j9oKcLCz0D9+/GjRi0EJdTa7akwuQmNMY2+8X79+FhYErEa5JQXOv7948aKF7UpOht8NeyjYVcA9urD9gxDWAjxjzAJXmgHX3lEIrnOtBE4vbRe4x64KTsLanjmDrQU+uy44G3NZCcGVt2/fWhhSfv782UiAplxTb2ARwcFYCLB7yKYUCwUdkMkJAODDZTsAulBMZF0NHn5+CAt6FngR9j8LzjlwfO4OVqtKWXCqFOaSOuBloGPod+7csbx7987Ci52/f/9uJOGSWOBs0jOh+jEL8HCjh9KQMHnxWTzioteyV6PwnLFV4GB7cI8e29ple5dUxa4FzgSTB64NLsbwLGwPzvBAWHFqxcrvMtGB9OzZMwu3ogfPQvdvx9IO3xa/pVCEXYTuj8U+ktrKu4zPnz+3fP361fLr1y+L/kvCyUUJeyWNYkcRcN5AUqqAs0VQF9y/H8hXPML3vz12EXgMOws8PA7YMXB6NuDWw6uCa2zz0GE8eIidBc5zP+YTwLklHz16lP57j50H7o/psauAh+isyn10LHZJmad4czYKzi9p7GN4YcxjvKWmZTIpA+4bnFdd+Jd2qoL36NEjBfdvHlQBz8MGfObMmZlR/Q84x6NKef36teXLly8WXQQDZ3Ly71Mr1NKcaNZLjTH0ImjCatWvWHnFl00fnTznwCt3wg7BPbpWmexUcqdSGJTFVtQBWwquHlwXnJPy2DSGiqBZcE6+CjjL+rrg/oGG2p8HrgAebo+EbdXdm2gSoowB3D+hV7gAeScajouA63M4IKtMj+xLOCoLJscnT56kWNTsrAUA540o1gHajw7B9TkKHSv89kYMGux/Epy0E1zQrPTY06BA4Dg8MMjDVlRI5IGzNcDGHNvUTKyU1LoDE30NgtqYA3CFmABZceo1tqzbsAw4e+7+bVQPztcyGLv1M/s6nIt/vzsEZ8igZ9cFD7FJHvY/A67qRz0YbMD52YOrh8fAw9KThx7Av3jxwsKdw9Ci3w1Xp1nYMXB/Lgovf/IVHB7gqLJjMk9UoDNJsTQFnAmJ/Wo1jvc5FLCF6PfCFUpJ/RxeQB/AFRYMwmGGD5/YhOjtBAc9hG4LOH9Hr2XiqALO7wBO/f3q1atccI8eW6azJcHnstYAmnWG368hDDcxbJIHzqsVCtva/k22RA2kB1CC8ayPW0ffkOCbyDzRzgLn31IB+K9jCAFsLSjKgIMegitZm1D/O3BVH/qZlSP4QLBUplyiIToXxm6tLn2qgLMSZY+GGp1JTQ0HPfzWGRWZf5knTBa4x1Z0DsxfDHf/AT3t4WDmuiQpAAAAAElFTkSuQmCC",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_a89a77ffe85cc350b6d5e694476e9cb5 = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAE6ElEQVR42u3b10vsQBTH8f3/3/VFsT0oCnaxd8XesPeuiBUbKiJ7OQd+4Tg7k53JJpvZvQpf8Hr3pnySTGJukqmpqcmm1fr6uraFhQUvm5+f52ZnZ7nh4WFuaGiIm5iY4La3t7N7e3u/+v7+5jI+gvuKHhUc2KmCh2H7Ck/YLS0tQc3NzVx3dzeHDVDS4D6hS+wwcEqHXTLgtMf4AF5dXa2to6MjJ6CXLPjU1JSX2CbwgYGBHOzUwF2wAb6zs8PoacCHYevAOzs7GZyik2vJgxcTPh82RcBqEpxaWlriSho8afgo0N6Bn5+fxw6eBHwh2GUDboMeB3yh2FRfXx+HZSpr8ELg48D2BpywiwmuIhQL2ztwV3QJXii6aSPEie0FuMT2CVwWF/YfuGOFYqcOvru7mxp4FPSenh6ObrmWDbgLugqe9F4uwefm5rg4wGtra5MHJ+y0wV3RdeBbW1scLYvNNAhXV6YY2OUGLjs7O/uVCTpxcImdBHiS6IWA50NPFdwW/Q/cEft/Aw9DzxQD2xdwW/SSATdhJwWeFHoc4Cb0TDGww8Bt0P/AHbF9A8+HHhe4Dj0yeHt7e1Ca4Ens5V6AS2BX7Hzg+dCTAA9Djwv88PCQ04JjJv39/b8yQVOnp6f/Nfj4+Hh2cnIyFDsHvLe3N0tJ8DBkie0L+MHBQazoXV1dHFkAfHl5mcN/ChM2wKmTkxOOTCQ2dXFxERSA4wlQusOVD/vm5iZ28DB0G/Co6CZwHOEu4DApa3Bgx42ugs/MzHDT09OcCk7zBvDR0dGvtOCUzVDy8PDgNTieQ4zrikVixwpuix0FnJ4itQE3obuCu6KHgUtsEzhOqLqd6/r6mgvAV1dXs9Tg4GA2bG8Htiv48fFx8Kx0uYED2gQObCtwiS6xXcAJW4JHRTeBq9gSfG1trSB0dV1U8I2NjZwwbwmNsG4B+OjoaM41OD4UBRzYKrgNug04xscwcGQLbjM8qr/4kIXa5eUlh2tyuW5GcPkhV3CJTdGh5IKOZUISDwGcPp8PHGGZdeC2J38VXIduAudffOiaUj0kdAi24Co2wF3QbbEBLtFN2Pv7+8Eyy1ywTeBYhqurK059NALvAEUGv7u7C7YiVlwHrYKr6EBSf1GQex7GzsXFxSAdOABpui7gLpe2RQP//PzkooB/fHz8ApfoJnCJXgg41ouw4wDHEKGCj42NceptkoaGBq6uri4oI8+4mDC21MvLCwdw+h7o9/f3HD6L5FhGXwROn5fgQFVPeDr4fODqOLmyshJkC047hi22t+D4ArhEjwI+MjLC0ckdgBSOFroOltgSHFc2JmyKzhc22Ki1tTUI723izmtFRQVXWVnJ0TOK9fX1XIbGGCw8Zg40Ffz19TXniuXr64t7f3/n3t7eOB24ek2KO2zqMARsDFvABrhEjxNch65ecdBnab4lC44VonAEyJ+ZwAEEVBqzcZLC+zVh4LRB5UbG9OjeiTxiTeAUhji8fQzwqqoqDvfAGxsbs21tbVxGHp63t7ccnRSpx8dHjqCfn5+5p6cnDhuDQCm8h4gN8PPzw9H39O8pTA8bQIWNGxwvNeHn6nBlApfoJQFOezodEfg7FRx7OuZn2gA6cISVxdAhwXHSUt8gcwHHn3F00LkO09vc3ORwtGBZmpqaOFwG0o6BEyquuv4Bbij46rTk9coAAAAASUVORK5CYII=",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_9be1b876fc220c2b22bc75a44a970ca5 = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAJCElEQVR42u2cZ29USRBF3///EUiIDyCMENE2OWeTc44mJwM2iKBZ3dKe2evr92BmbWu9mkG6kmde6O7T3VXV1T00hw4d6u3bt6+0c+fO0q5du0p79uxZor1795ZmZmZKt27dKt25c6d07969/neXLl1apCtXrpQuX75cOnfuXOns2bO9kydPlk6cOFG6ePHiIl29enWRVA569uxZ6c2bN4v08uXL0uPHj0vXr1/v1+XMmTMl2vHgwYPS8+fPS0+fPi09efKk//fdu3dLN2/eLGWdeeeFCxd6c3NzpZ8/f5b41xw9erR3+PDh0vT0dIkOkPgO0OogSZAkGn3//v2SKgRwFSydP3++BDwazTuOHTvWf+/BgwdLCV6wJDrWGw6I2dnZEtCAyDPqZMpE1JF73717V6LT9L4x8BUETt2oE6Dfvn1bevjw4RLQzLLTp0+Xjh8/XqKuN27c6Hfc9+/fSwsLC6WGXpJ48MiRI/2GY2a2b99eYgpSURrjgOlpTAb3olOnTpVUDqZq//79JerCvUDDLLi4huh0OlwNlwQJc0bdBgEu2MMCp8yRAA6YtPcSQHO28UyO7NevX/cePXpUSp8DaNoBN12jk5hlX79+LWFiGj2YwGm8BOjJyckSI4NKeKP4jl6+fft26dq1ayWuU55XlsoDho6U05JwgC9evOh/R+PGwFcQOKAB8urVqyWA8SN8xt5jRgSWejND0yfovRJ+QKJjc6ZS13SaYomD1cCQAE8kVU5TTkvS30xxAGNa6FUKZ4TIhCT8dFbYPqRKY3YQjaGzqKgDB8owwKnnagOHz5oErhGcjWA0JnBNbcnjagDzPpwkn5lR1E8CeFdcj/kQaOBTds5U3gVwPY9fATTP0nlcbxTq4Ri3bNnS27p1a4nvdu/eXWLq0znuGGkUoyfDP2BS4THwv2EvB7jH1Tl9gUrjPLTK1SffU0Eq7qYEAIBG+U6J7zAhHipKgP748WNJIdznz59LdDr3Ykqw2QBXJ6QJxIRglrjebNu2rbdjx47+cp4FDs5S1yWFcBKF+mjMhczvgGcsOwjwtN2rBZyYeRjgwM4ox6MqH/G/BT41NdVfaQKcKelRB0AzXMrlsyDRwHRwXWB8IcOURukbgNyWt2FQUH+eoTMFhVEOLOpCedkuDQg6ZWDgcoha1uMcCQcFe9SB52zkXZQv4AB+//79IrHS5DpqyGEAmigF0AcOHCgpgpG8sZgMACPAc50KC1quFtMsZAd4riSFaQEMn8m1SGnW0qHTsRKjcc0BpyLDAAfoMMAzOTUIcE9uJXCHhnQPoBEOj1VjW5IME8Gs+PDhQ6scfJkUgJIdzLw3DoIR4jYNUwHgzKEgNxVAIrTzXIwnnPQ3DcQMAQDlEt+X+SSyMvwbA18mcM9XS1zLWD2jFspR2WlKcIR8TyyNYxRIQGf+HJOSwOUnKixkRUm8DXCcJBDdWWH/clGRy3QP9XiWhuZyHBDc52YGiIAHLiDGwA14ghwEuDsyiXuw+0x5z1fkQiNnTkY4mmWYHUADCZhofn6+9O3bt/7fAMa0pHiX4Hs+pQF05k5wNoxSNyV4eu4BeOZb3HZmg9MBesw7Br4CwN02+6ZzRhp6P9M/6wDMXC5LTHvAp0nJvItvGAPnx48fJVacwCWn/evXr/49mVFMs+TZRzpW7W1YwmNKAA40XzTkHmWaEq5zf5spGVXgOP1lAc+Vpe8P+sLFG5rRCB3StjuDuAenSSPVmHR0mB2A8ywmS2VzDQdH7oQddqACWvBlTtykYDrSweK8VU9GtkdZS4ADMdOZEkAZeaMO3KMR3+WX3Dcw2tURS4ADk9HsSuCYEtKUdJLnMnLBkTkOpA7lfXQg5WJ+WNQwimQ6aEwejUiT4llCOp1nP336VAIanYZjdPg40gz/BgaeTnMM/J9MHzE0RxykvId38BkzorIzkyiT07CibItOcg8PcU/uYme87jG7b8n57pArox2+z8NGGn2Yl0GPTuh5nsGxAolR2wbcYQ8KnPqOLHDf+vKNXSKNzKUA1c8p5mYw8kiJM4m5XvD4vCGHwhaa79bwN9fy+BkwPW+SKUw/quByU0U5uc/pCSbJl+2Ek5lvobH5eU0DzxG3loFj392uSzTaDxJ1wfOjEal8H8KEZJTUtsvjm9INm8JtUzudWJoUvz8XPYRheSjHnSnXclcGR5fH3DztOgzwPEy0WsC7ttXWLPCMbNI8uInI5FdmFBnZDrQrG5imRvd0rVyRZw7JmXTJ8+NN7jv6kj13qRO8nzzKLSxCrzx56lt0ubwfOeCZI1kJ4AnGnR3X8vgD92I+MvJw05L35CaFH6fwvLmbHzoHEyCHmZ2RUVDmYSQObKLMOkoN09mX6gBmHxOzk8A5HidlhvB3wPNM4GoBz/MrgwAnQlkJ4CyolgU8j+vmee623/w4cE333EmhMb4sdoieN8+pndGIA89ceK5UM/vov57o2gXiXt6hzgHmly9fSuRbWoEDU/JjxFKmXdvOR48qcD478BzhpAa4Lpu/CHgCzvPeyHMgGbH4XqED98Z2pTZpXDpIgc73/Qm4/5oOM4eJwZSkiXEzkyFjdtK/Bu4mY7WBU+lhgDOyhwHOtUGBe3STTjeTYXz2qCfjbZ7NM5di3bCKlDAVvq0mZSU8vs6fY+QvGDz/QQO7FhGZ+VOjKDOB5ywAuMrKuuWuU5uDHQNfJnBPUeQKNlep/neX78kz6O5rsAgkAjnqjUh9F3A/dJmhHVB5sReSFckVZptNTDvLu7g3HZQvqSkHWLnEz5+WjIEPCTxzGDQ2Z5bvr/pslfKHABnp+JZf1++F6BzMkd5D5+av2DZv3lzasGFDadOmTaWNGzeWCrjnnVlaUxiA2pJHea4kR4inKP/0Q9Y8UeXiGiPr/wqcrcz/HHibo3IYeWQ4TwywONMZd/+5o5Q/8spy2kxI11lHB44TZg1ClLd+/frSunXrShMTEyX9hKcP3J1jeuL00K4cEblY8VHbBbwrMhgDHxI4I74tB9L23rYjb7m5QeM488iWIGfXHTgnEDxC8qyjm4qsU55Tp4NlMnCODt//KxEsBHXz//YE8/MX7G+51/YmauYAAAAASUVORK5CYII=",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_d1707f042f199bc480911a410925d0ed = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAABtUlEQVR42u3ZzYqDQAzA8bz/M/UderIexFIQxC+wFylU6C4OzKJbK3ZM4mTM4X+SufyETGDgdDr9SK4sS7TqunYqy7LVgWK7Y38DLR58T2wXaAV3wN4CLRqcGxsDWiw4JzYmtEhwF9iqqmbjhhYHjon9CZwS+hDgc6OjaRrT/X6fxIEtAhx7Vs+Bc2F7D459Ob5er7c4sb0Gx9pGLpfLJAUnxO77/g38er2a0jQ1HR58C+5cCo6M/QnaFkXRpDiOTUmSmG6321+HAaeAHhq2kTH2+XxeBOdAB2nYz+dzscfjYbLr3xjbtgROjQ4SsLuuW13btiZ71s7sMfjQEjglOviK/Q2yC7jtEzgVOvgE7YpMBU6BDntiYwBvBedGB25sCuSt4JzoLODUyBjgXOgQCjQHOAY6hAKNAc6BDiFhY4BTo0Mo0NzgrugQEjYWOCU6hAI9B57nucm++Nhv34aJDtKhLepSW8HX/BAycN9x9wBfyhlcCq5v4P8DH6AxcUWDS8QVCS4ZVxR4CLhiwEPB9R48NFxR4CHgeg8eKrKC71hRFCYFV3AFV3AFV3AFV3AFV3BPwH8B7rUxKEpomD0AAAAASUVORK5CYII=",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_047cf7a4389ddf495e7be954fc800a97 = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAH10lEQVR42u2baU8UTRSF+///FxNRIwoouIviBuK+4K6gImBCQoK8uZU8k9OHut3AyBuQ+XASZuju6X5u1d2qurlz587OrVu3im7evNnSjRs3dq5du1Z09erVlubm5qq6e/fuQNevX68qruvif/x23Fdodna26N69e0X8zv3793cWFhaK5ufni/h87ty5oqdPn+7Sw4cPi+L8EJ9fvXpV9Pbt25Y+fvw40Js3b4qePHlS9Pjx4yJ+d3Fxsej58+c729vbRX/+/Cnic1N7eKBeuXJloJmZmZZ4cEAAWQ3n505PTxddvny5pYsXLw506dKlIgyNuCYGCGM8ePCgCGgIA7ghHj16NDiH+0ac6yxGwP8ScI6N70LA49mBquJ6XAPxPcfFzNjc3Cza2toq4nNTcyUBZ2pqqiWH5FM9AKiBQgBGgAdqXHdycrJoYmKiiN8DNIMAt6dGzoADUcX/Dhs4ruhEAAee+3s9h2Px6c+ePSsCKu5jaWlp8D9cBWDx/4gZFK4Et/Pt27eilZWVIj6nwAGBAA08RkIWVFXuHriWXh/Q7lK4J4BrUOaBfTSPgA8B/Pbt26mACQgHrsYh6+EzkICKNEiSuXiQdOAY5PXr1wPgnz9/Lvry5UsRnxtN/TTA+VQHkgdPFa6kD3gc22dAjH9cgAfsIwtc/Tvn+rHk5biSmoDnWYnGFYUewkjAe/HiRRFwP3z4MAD+8uXLItwP4DGsBk1Ao/fv31fVxEPy8DGlAXzhwoWWmPIO2QOjjtYR8A7gAD0o8K6R7xCBrnIjeeVZE67Fixjg6ncABhqg3717t0tkRvhmr2i5FoaIYBmGCuE6EKBxP03AjYLj/PnzVY2PjxcBXKF0Qc6A6yg/TOD+3WEBD9jaAnDgq6urLXUCj//hy4FH/h3TMqskMxej2Uw2Q/rckc4UdztahWalu6eQ7i5C/Db35v0iwNd6Lr3APUh6sDypwNUfe8+EcyL7+Pr1a0tejO0CDjSt/moVJsfplAaAG8uvwbVrVahXsFnmo/Lejw6CTG4AfCrwNHB6ADyWwLOyvwt4DfYwwDXQoXADwMINOPDw0yFy9TAO+bXLZxAVZlw/fmvgUnhgwGddu5oyI6lbyYD779V6L26ArNyvdQwZlZ5ZjIAfELiObOAqYKS9alX0TAAONAe+vLy8Sxz76dOnln79+tUS1x7k4TyoT2kebi85tBtrL8Azt6HXys4lHXRX8s8Bz3Lovw3cR7hnJYDW6tFBA5jMAxdCzyNSOSCsra0Vra+vt7SxsVHE/3/8+DEQXUBAYwiOBXgc06o0PZPISmyVp25ZWriXvot3G2vnj4AfAHjWus0yka6spNZJzBYjPEjiJhTS79+/WwK0K9K6nz9/tuQuBOAo3Em4rVAYufECp7ZynwHvK3y6euQnDXhkNhEzhgLubkLz6WyBIpsFajjPzTMXpZmSB1DvoajLoRwnDcSFOCyVFzAOWl1IiF6NqgD3lM7h1lZzMr980oDXwOOrdflOZ9gu4Exbna7Zph1XzWc78K68G8BZa6ArRvh9Y4Da6jwjjn63uwkfxQ67a6T3AvebP+nANRPRpTt6L74jzKtU+jFZTdD4NgRuWNfs+K5r9SWkYLOq1N1FrX/T18+p9dF5Dl+7JD3UUcZ3QKCIUdgAx1D7Aa6wTyTwCFYBgurSl+RU2o+JWADYTKR8KDIgn1Xa12l8Q4/vJtJ9GN6B4+GQVqPZxs1s65wWQ5nr0MDt18Xo3CvPAeQjDdxv4igCBzC/zT1RWRKw8NMBgyqzBhdlYD2ecR9cH0VsqGVGVL+NLxup3/LVacBjnK51Rg9egNH8PlurzPot6pr2A5zRd9jACcbHBjjHcG5tX3nWC+dcHlJ7KL5Vgf91Aef6/DZG94CIYbUD6UG5FTRrPjvbaep7q2s+NDuHz7rvz7c1jID/BeA+a7jRmHpMcQD4wyjsbLubzxLfYBkP6cABkwHv6vVkKZ8uQmfAS9DUqO7B0Xct+a7UAOzA/UZqwPF1ewXetb+wCzj3sB/gof0C933kvttWGewbuKY5Kh7OV0BUrO+56OD5FuIu0Blw1UGAZ30jql0Hrq+aZBuHWsBr++5837W+sjEM8GiF/p/Aa+rLq92lOHBPAkJZkuEeIc4rwLN3VTyp100ztQ00uqHdt+vqeqAD95eWkAfRWjDtextuGOBdxxwYuAa6owi8Nqr7gHeNcN2h1SV3p17kKXQ3Ap81BaboKyPcXyxiStOYAQwQwz1oG1LFLlJA6zW80YOr4RwHHqPDC46+1af4uwu29ou8VhgBHwK4jvQMtDbe+oDXagVt2oV49dG3evN97EJuNCh6K5Plou/fvxepOwAeW3wxAO4B4KR+9IlV/A8BPoTRvQYAYtdyn/vxvQB36McaeA12F3B950aBxz0ChBGv7w0hQHv25Rs5UUDzYiV+S9NFDYCA5h5osp09e7ZobGys6MyZM0XxGjp/N1GAsMbn++R8yclXREIeADESCiMB0oOwG4vj4m+Ae857GMA1h/bG0z8J3MGHqEa9UtVWcLYbLFte87eSURdwf/s4DIQh9TXL0OnTp4tOnTrVEgYINRGceGAgAbMGOltwdeGG4noYEJi11zxGwIcAHt+rv2eEuw9nZtVGPD1mb4YBvPZmXAa8TzXg3guqvcvvBQ7VKEGytlf+P1IB94zuRhRRAAAAAElFTkSuQmCC",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_35a70477cdb1aec56ce558f46abbb54b = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAIBElEQVR42u2baU8UTRSF+///CmL4aNQQIRoRRXEXV8R93xF3UNzgzb3JMzlzrGp6XkAxMyYnGafpquqnqs69VdXTzM/Pb5w9ezZ15syZPp06dWpjbm6uT6dPn05duHAhFfer4jv+Ju5XUW7pOm24ePFin6iH6zdu3EhFW6anp1OHDx9O3b17t0+zs7OpmZmZnviOtpw7dy5FPZcuXUpRj+vmzZsb+/fvb9Xjx483vn//nvr582efmsuXL/8GAsVDOTQHXgLEQwAJ0Wk8tIpr3Iv4/uTJk6mFhYVUPDwQa23hXgWOeI54/tCVK1f6pPWMgG8DcGYHgK9fv55aXFxMBdjQkydPep9RCbDr48ePffrw4UMqgQPE4SpgV+0hQw7NO+Do0aOpsIFDhw6lsAXK984CyO3bt3u6c+dOChuotePYsWM9wEeOHEl1AR6wR8BbgBNH3CZi9jpoxD0OtU3Up1LYr1+/7unZs2cpOq+JBgFEAyYP7oHUg00JOJ1REw87MTFR1fHjx1NYCWVrUAQAYP8m8Pv37/fpnwV+7dq1FN6t/oqwgatXr6a8A6KdDhqdOHEi1QWwZ2QqB/7o0aOUW01aigcoDVIA9imO2hpR6gw65ODBg32amppKUS8PqcA9oHUBzrPtJHCHveuAxzUA491MeerDszUruXXrVqoGHNC0TWcoYB24lxWx4d69eymg0ZGlNUco/taBP3jwIAXonqVEIyKKayRXETxrgbDkZ12m4NADr025rQCvjXBGeYh6yFyYYYAGsoJmEcLf0DnYj87K8+fPp3zFWrKBUEB5+vRpn+hQyvKVeVyrAX/+/HnqxYsXqSYCUwD3hQGjHuBto3U3AgfKIMAZhYMAj+/Dlinj8+fPqbdv36ZevnyZImtpBR6jn8J96RtTDLC1bIQGlsQKk6nu+y10SIClTqZ4TVqvrxNIFek8bIMUMz4z/bnmVlXK1ii3M3B/cNewAvdnVQ9X4MvLy6k3b970CcDv37/vU8MUruXbAKJSFQ3wzMUzm5L3A8RB83/apdBqqRaivrgHe8GGWJUigGtwY5SWBtdQAPdR2gU495SAx0h2v9aARoDzncJS/h9t5b42sZcSnRDxoCHfrsEBeGkbE0up2cZmK84QluL77iHfI2dUPnz4sKgInFiGj2isY7uAh4YSuO6tlMBiGe73ZCAx7bEFdvSYMbV8n/+/evWqp1+/fqUImr5bGHaSK032LNT/SioB57tauqj3uB1xzYNmabeQDnT/9ROePw08AO8ocA8mXYB74OkC3GcIM8ltQsWIr/k8FsJUjykOYIdDSoc9ARqobVpbW0t5maGwvB5wT8dcDluB1xRBSxcuKkZiKb0bAd8G4D79fVdNBTy8G+A8eCm1q2UuBELSNexDLQQbWFlZKWoQ4JxdApn0kJn9m6WUwG4mtwvdUnXQ/yLwEP++fPmS+ivAKcR39jSlAgryYFZL8ULMFOrhaK20D1IDjpUAOazi3bt3KRYjQP306VOK60tLSyn9B3Cs49u3byk6ryRdHDW+hGd0qkqLHtWwAg/Qq6urKfZMVJTfClyh1XbjXLV3OEK6V1ESwErXKJ/6Sce807TjAE3ax5kilhIA3EoAo4e/mi4ymnVEo4GB+9J+2IGH7379+jWlL/DwYg/gGemURScpcGYXbUzgvrKsLWNDZBz6qkLI0yeVL60dqgZSv0ZDa4E35KAJkroKDMXDMhp9tAKPQDcCvkXgATG8GUiABiIi01hfX09p1sG9+Lzn2Ph9AKYtZF08Rx6x+amNvkPnoGtninoeOEj6t9nuH+BLCx3K99MZRheLF4DsauC1kb2bgDO7GNF+sgIYfFpTuR8/fhRVs48QdlQDTozQ9jMIaGMCJyChEvDaWaK+nANYT+18s4ggVhINQ4wMtzBVF+CeO3cBrrC7AvcBs+uBYwMIDy7l6Mw2D44AJ5cupXKMYEC7tcTfOHACaQ141EWne5s0iDa1l87bgHsnBXRA+CilMl14sPjgs18bAd8CcIKJ576q2pIY8OzwacrlD0U92Aerx/gMLGyhFjz5v8NWUS65uw4a2kAbsRYGR+6llBY6m8mBhyi8BmKrwH0V1wU4nwcB3ga7BpwZ6sBL7+MMDFwzBM0aGIV+Sg08Tb1IuzzlAogHKH0oFxbCw1Jv3F/b+aNcXwiV5G3TNvnIhgGsOgPXNM/BDjtwTwMjEfA0F36lk7NGAQOPjMJPoPle96wRf1OzBwIUsEvAa2AIgCV5ByM/T9RNJO8cLCcshXYC3O1vy8A1KO40cIXdFXgb7Brw0gHuZsD912aAp9MBrK9U+AKNmc+2BpDZIIyf1jS6N8LqkEJpGJW5V6m4RgNdGrxcTFcfTXEfUDYDrdc83SSI0UZsiOcaAd8G4NEmX1jRRn+uqMeBI9pImktZMdNZCQOY7WTOFngHPhTAx8bGNprSqtAXDyqHjxhVtfNBTcdq8rQtPnu9Xj+jtzQIdhPwgP3HgDvILsD5rPXqaQz5t7+mpnLgfpAc5eLzBEm3FL6nHt0NhBsvLwH3wIEDqX379qX27NmzMT4+nmp0f4OG+ql2F9HwNrg1S0GltMw3+EfAtwi89MKM/nKXGeR7zINIAStk9sk1J6/FE8qIDTpfj3B4A2h++sjvl+Lz3r17U41Gax8BXQC72oB7IK2VMQI+IHD1bZXvkbjUp90y/g9onqeU/fh+CvK9En3BnvSvdtTI75T01x38VmpycjL1H+2Zxlqr6oJTAAAAAElFTkSuQmCC",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_aa3f48d500ceccff387b5b78900cb83f = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAJNUlEQVR42u2c127cShBE+f/fYMNPhh+c4JxzzjnnnHMO8BrVwFmU6pIrUpINXckPBYi75ITTM909M1w1V65cGZ05c6Z05MiR0sGDB0sHDhwYbd++vbRz587Srl27SlyfOHGidPny5bHOnj1bOn78eGnPnj0lyufzkydPjsU9lM916sKFC53y8rp0+PDh0c2bN0fPnz+fVrTTRX/R3r17xzp06FDp9OnTo1u3bpUo6+XLl6VGwLs6LW3ZsqW0devWEvfQADp76dKlkoAfPXq0tH///pIDE/A2EGnQPwUb4H2gtwGfBH3eAaeMLtiTgM817LmEDmj6q/qePHlSevfuXen169el5uLFi2OXQqP00ObNm0vr168vbdy4sbRv377SuXPnSlevXi3JcLgUd0nS7t27SwCn7Dbg3LNjx47Sn4Q9V9AdOHXOS+AoYTtw2jZT4PRpEvDpoD9+/LgTukRskluhXc+ePSvhSiirAZzkjaTDAMe10Hj5KYlnFSSxeAZFPu8CLjlsCbc0CfhMYQvAEOgCPgn6ggLOLKF+dO3atRmPbAHQrOwLfTrgSAxu375dwqVQxtOnT0uNg8ZKEukg4Pn81KlTpWPHjpXkuzKAyr1IGMNTp77AaVMbcMGeBNxhdwEfAh3gun8S7HkJfPny5aVJwB2eYG/atGksYHcB7wsb4H2gO3CJmMQsPH/+/OjevXulV69eld6+fVvKMhoFN1Kvbdu2jTZs2FBas2ZNiY4SLKmkLVpjDKzNvUBuUwJPgA5bUm7bF3gbSAfeF/ofAS7YswHu/laZj0T2wSxx0Hfu3JmiNuAsHlIAJ47Qjj4ZSRfwNuhApD4GElkY0iwG6Pv370sAJ3hiiEYBce3atWPQiJEHLFwMHfJsIQNcF3ApQbumA90XuO5pg602JvAXL150Qu8DXLDv3r07BvvmzZuSypWoZ5yHTwLuoJlGQLlx48Z/8m5WiYCX3OdOgi0xSpA64mI1i6E1lSVcmhtsCPCETnsJgCxsGFgYWn2nbb2BC4r7Yfx5juyZAFe5Q4ALZkJ38LMBDqyEDhhlEA5c9RL8SRCohxmsNtHGzLcRrmTsUnJ5ClhAA5HGAlmdoMO5onTgEh2RkdwALjUOoAn84cOHJQwDCEZaGwiHLSP0AY4Lod55A5yODAEu2F3A79+/P25gG/QhwJkNDjuBA5kpz34HUKgXA2Tb5WbIq3NFybXXISkn13MNO1xkIUAj4AGPzjkoLI8bImjRMb5Xh0kZXQKNfFQAHXj4UjrO98DkWkbhGYzSBlx1LDrgDjuBA70PcJXP94LDyAY49/JsrvwEOFO5dAMZaBX88h6eJWjms7isJpfUgGZlCDxSIj/ZwQiAZmpz7Rvz3AvEBK4R5cDdZfFsCoMuWOAOuy9wPwXpCzyzEsCyK+fuQfIAmdOfa/w/AORC+Bs4gP/06VMJVwPEz58/j0Waxz1fv34t/fz5s/T9+/eSysGtPHr0aNSwkmRjKE90GGVKw4jOKH1pnokq2OKyKI+OtwGX5Ov+AZ8D4HwH2AcPHpQITID3bAQ3Rt6t59Vo5O6B8tLf5xqA+jX9Afzx48cSkBBwgSl4gMYlffjwofTjx4/Sr1+/St++fSvpXu9zw8oS4ADKs0oJAIjDAUZgBs25Bi79X4Grj+rbrIB3bYnqM9+SlZjyLIIoEyBaHZJvE/AwHNu8169fL9FZz9XdYBLgiTfsw6jj7FXjShIshsBdfPnyZfwdn1GGvpMwEtcqlzYxqAq49lEkTnZIC3ElgpfZASN7sQPXLOFZQOPvNZuok9y8DiASOHkynRRgPgM01+mzyUjkXsjfAc6UBjwgErwEaAzIPRnc5EIISN4pCQPkBpiA45J4JoMlKR5lyy0AFODcMxj4qlWrRi4HzoheTMAFk2vg6XOAAhrhw2kbscRnos/IZt26dSNEoJGUheShMNfAzDeQ3Cj4eVwTMyb30GUc3FiucjEgriQDrpRwgAhUXIyeyy3TXIYzer1MRv0Q4On6WoFz00IErs4qm8nVYO6DABw3ob8BmvVwLwYhY/IsC2Mzu/RZw8kOnfJTlVxa535IvrwpkLkaxR0AMV8SEnSH7ate3BFwASToAM1NJK5zCT+vgeNL/xZwlNsMk4AzE30fReJeZgOG8T3rTO1wB8BMw0i5n4KLSVej73zVLAFebBtcwKRDgnQLfM6zbSf3fu7nwMn3ue4DHNACmUv2PsDzkKAPcIfdFzjfzRvg/jJPX+AedNT4tlWmlJtidJYZIWE4Vpa5Jez7I7mD6CdDEq6EsvRMmxvza7WnYcqTrulvQHou7gse9lL4XlAZjbkIktzF4N/8jdkEzohACxY4f/8N4G1K0Ck1HndAJ/gOF8MqFVfjq1GmfAbAtqDJPXlvGgkjYmAJ4+R+Ue2HA5kXMn2PBOUbVVxrVHJSxNlnvjAkySCzBU4HhwDPw+g+wMlQhgB32A6cWUcbp4zwvsD9fXLfIWQPJYUfZ1YocMwEOJ0AHoDzEITVr+rK18+6AiHfK/0jCOIyMADXPINrk/vKLGoicOD55hRBkkCa++CLGTjtIFZoJhH/MEJuN9AOzfCGLMPl755I5Oh+suKHAqgNOKtFz1pmA5yFGW3ygeDA8eESIBiJ1AcIdym5lPe98DkBrlE6F8AZ9fl6whDgk6A77C7ggqk2MTP7AvfRnruArDSJGZTpx4T5enbXYboyskZQ0q0A0/eq/fUxOivxbAIHgp/WzGaUO2zfS0k5cO7P89PcxliwwN39DAHeF7YD95OqNABtpD+ChhsgK8GVkGfni52KH8wQQAOYH6CtXLlyinTuUL/TzFMYXEounwU/07t8aTLfB8lfOaemA+6w24yR9frh9qIF3gW7Dfok4G2iHsUepm7+uJbYRBtxDwJMhkJ+z+Evh8J8zrMCSx9hQT2AXbJkyRStWLFi/AOHxk9FcrRhQdJESRX4ygljJHRBmA50G/ghsP8BnyVwCdgZG7rkv4Jgq5nOsepF9JMVoTIXYhH1Aphgma9FaA2QP7PE0ABfunRpadmyZSW1hU27RoByOpM/AlMjCeAOuwu45IFsCPC+oBc8cH/NYDrgbdlDX+C8VDRUDnz16tWlLuC+oQTw3CvJF5R8dZubXfkrC97tkfIfH/wGz3awyBSztDMAAAAASUVORK5CYII=",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_6610ffc3843fc286b5f922d0ff999b64 = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAIcElEQVR42u2cSU+UTRSF3///G1wYTFgYiIkrF0ZBREUcGJzAWVDAEUVAROwv9yZP5/SxqoevNSHpXpyku3mHqqeq7lSlzfPnz1uPHj1K3blzJ7W0tNTW4uJi6tatW6nbt2+n+P3evXuptbW1th48eJBaXl7uEL+jlZWV1P3799vye7jm8ePHqSdPnrTFb7T/xYsXqdevX6devnyZWl9fT2nbeC7fuYZnbGxspN6+fdva3t5O8VswC3GP8grFe7j23bt3qQ8fPqSaN2/etP/48OHDFDcGyPn5+Q4BHiDcw8tDgOAankcnV1dXU/xdB3RhYSHFAHAvUBnUAO4DyPufPXuW4rsOEu2lDQiIAEJbW1tj4H8TOM/lWgDFzGZ2w+fp06cp+uGgdTXu7u6mvn37lvr69WuqiSXHiPHyAH337t3UjRs3UgDne2lZ+3KlIYClkzw7TBOgEaaKZ/EeQEeHMRUA0EHQa1G0jzZ5W7oB95ndD3B4jARwvvtK6wace1CseFY99tz9Sm3l6irHhHz58iX1+fPnVBOgfQnGkr5582YKwIglX2oAL2bZMoAuzEU4aOD7M9wRAiJsKmLm0e4x8L8AXJ1XKD4DmpkOGFYq/WGVBFgHznfuffXqVYoIJ7S5uZnyaIjnApzBi/sxQ+/fv08Bug1cZ0SIMA8Tcv369RS/o1IYR0N0ZqnUcfnz3PEBU4HzeRDgwPjXwLn/VAIPOQg6A7RSPIx6gUbxPFYbbeA7wpRoGIc5+PTpUwoz4z6KNkd/aAv34CT/MClxMfDCZFy9ejU1NzeXmp2dTTl4Be7hni7lEJ0bAxfg2Oj/C9wTCZWbi+ggAHxQuAaIKDoMaB8UVHKWunpVmADM1M7Ozh/Cb2BaPH9g8sRAYUL29/dTBwcHKcJCwDcBWcE66MuXL6cYCLW5jHA32A6cGTUIcGbYaQQObI9KXO04vBtwZnWIGgqQNf6lAdRiiDzc1GgW6LG/L32tf/g9Hoc76HgObXPQ3It5UgeJw0N+jSdU0SZmdt/AA1DMWMI9wkGFPYrAY2bzXk++dOV5VAJ4j06w7Q0mAuAedwMRmw0gDZN4hlcSAa4NphPAUkeqUtMCFBybF5GQhnaYIcyND5oPVgyCJlWnCrjHpf0A12XYL3A6Nghwri0BV3vPdQH248ePKU2uQv4M7omVwbXYewcOYBS/xSpqHHAto9RlWUsAuMcrfuostaCksbmbgICHs3JpiBjSvzEYwHJoY+BDAtdSql4DLJ/x/B2oCsWjE6/Z8J64hnswGQyamxJ9TrSj8YwSc+COr1Tv1rq5ymvXWhzyQlMtaYkZVQPtmwRj4ALcq3H9ACeS0LpGyGNnNQE+GDXg6uR82QMN8G4eNMzzevfe3l6K74iIh3Y2xN2A93o0cAMeIBkM7uEZ1LR1Q8BjV/f02Ecgq8bAhwCupkKrcy7dK/TYGfl+ocL37NDNAVAVtMPz+gd/V1NBIvP9+/eUp/Ionq+xfeOgCQM1EQH+zMxMhzwpwsSUqoOjDpz2DAUcmL6Zqw4TAQuwJbOBmfE43zNMri+B7uZIySaBjskADvBKAjrChPz8+bNDXB+DwszW9zfXrl1rhXB4zFI1KYDlGlL1MfC99szmOlZNrCqPiNKkOHDfJtPYuOYAvf7tpdiQp9buNOOzPw/5fiMd0TMjLhyf7xLFIGAWMB0Ap7TqCqB8BvTQwDElCtyTkVEADsCafOa73Sfj1GMWHU6TzJLIQ52Vn/fQiMXPfWACvNrnCY1nnPp8F9cyKBr60QnMAmaiJnV+LhxgL9j9AtdDRCMJHNgeQTj4kjk5PDzsEL/7AR/NSj0n0Oy3IaTDpChsOl6rWTsgNT81E+IhX1zjAwhonGQpafGkpARZdaqB+8w7jcBrVToHrjsu7vCA5dmimgxMiWeS3Mtg8X7NDTzbzVoK8HTJO0hAuzNTU+LhXy1NLzWiW+m0VAIdFDhw/jVwL6qdeuBeFQSq1i8ctFftfC9RnaIDr9U/FHQtHOSZuqOjzjFEP7Qu33jKXUo8/DxGSQwKD66VUPU7nz0jGxngtUxvGOB+9EDNgpsDr0fTAV36fHb57ribiRLg2rPimpqz9HAQExb98HbzHSaRYTdal/YzI+4sS4c3NZ0PMUtL5zyGAd4LUgl4rxldgz0ocPrh7fYTYjGxBwautRM9voYNxy47aKBqHcJ3tL0eTUe0Kse9HtJ5PUSdn8Nyh8d3DQN//PjRIUxMt90hTEjtsGkC73bo0v+hkRerRhG4bzgHWC/IwctPJgSvRmeyx8FefNETSLU4m9+9Tg08hebLGCAMjp7xqJmBmq0uhXncozvpeoYkgAL4+Pg49evXr9TR0VFqaOD671z+NXCfpf0AB8YgwGtxdQ14zPJeoHkW/SmdafRTDF7Ui92wRmsjXtlzB8gs1di5lkH6AXU6FnIH50uf63RXpleSoim3m4ya6dK0fwx8SOClA5Ve6dOB4XkAR/zuTlL7jsmgPMJJ40uXLqWmp6dT58+fbzVas/YddEyIxpXElG5K/HwfHtuLSLofWCv8A0T3FX2QfFAAEu2snWAdA+8C3Ddy43OsEtpS2kqj/Z50MSheAWTvUUU0ginhHt2mgw9JHyeOz549mzpz5kzq3LlzqampqTb0RvcZ/R8T+T6h7pzUCk5+7q+UyHgqXdpTrDm4MfAhgZd2XqLjtWNntI02szq1uMZq85o5gwlMdYY+GF5Lof/xTkwI/zrkypUrKaBOTEykJicnUxcuXGhdvHgx1cQDgOigNXOqVQFLO+l+KspnI863ttU1Bj4gcI3dFbTv9fkuTQyCVw6JmPw/JfAMV4ETlrmjBebJyUnr9+/fKZyjg/YaTQw+faVtXkfibE8ADun5eszOf0iadF2eUGWUAAAAAElFTkSuQmCC",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var image_overlay_ba0f0d14f5a326f6d0321ceb4226f921 = L.imageOverlay(
"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFwAAAAwCAYAAACc2R15AAAKXklEQVR42u2bZY8bSxBF5w/nSxQpicLMzMzMzMzMzMycSIkCfrolnVG5dsae8Xr3ORDpSll7oPt0961qcHLkyJHKrl27TLt37zYdOHDApO/2799v2rt3r2nPnj2mQ4cOmY4ePWo6ceJEqm3btpk2btxoWrVqlWnhwoWm6dOnm0aMGGEaNGhQpXv37lXq1atXlfr27WsaOnSoafjw4ZXZs2eb5s+fb5o3b56Jz2fOnGmaO3euaenSpalWr15tWrt2bZXWrFljWrdunWnz5s2VLVu2mNavX2/iGbyX98yZM8eke2AIk6tXr5qS48ePV9CxY8dM/C2Q+/btM+3cudOUB5p7Dh48WNm+fbtp06ZNpsWLF5tmzZplArjUu3dvU8+ePaukhvDAaRzuUwPGigMWAW3r1q0mdQKgAXzlypVV4vMNGzaYVI8/GvjgwYNNAO4o4CtWrDAtWbLERBkBD/AdO3akIzaO1GnTppmmTp1aJV0Lt7Nnz5ouXbpkSk6dOpV+eOHCBdP58+fTFsJKsJ3Dhw+bdJ+XYMd7AL9o0SLTuHHjTLIDJNh9+vRJFYGPGjXKRCUFiR5Hb4wWQk8DtGyS8nNvEeCCXRY417YscKlZwHkfoJHKTc9esGCBidGG8H/ACzj30AjUg/dNnDixSvqMUQUvmBAfk3PnzlUuX77cRhEewlKwEAVWSY0gO5G4l5dTyYEDB5oGDBhQpWHDhqWaMmVKlUaPHm2iRwqAPFICCPoHvABwgOYBJ2tgqAog2Q/f0eMpI6ABJhsAZLQB3jNp0qRUuo77fKygUbhu7NixJkaj4sayZctMsfEpY3L69OnKlStXTAKtvyXAAhGwBEsPGnEPwAECcKV/Uv/+/atgS9jNjBkzTPiwB07hywAHWhngvqGKAidQtyxwyVdy8uTJJoY4ebG3BYIYhadSzA1iRqX7ly9fbgIWWQ6gx48fbyJWSJQXYSETJkxI5YGrTjQkHYXGggEZW6JJjQeG10QriZZCQ3jgfMbDeSlQhwwZYvrdgXtLpBOVBg7IRoGrt3Mvw5dh26NHD5OfSQo43idRGSzFB0mJPNpbC38DOA5jZTDYASAQMGlwH7gjcKCiLOAxCDM3oLFgkqinqFecPHnSxOQHkECk93CdbOfMmTMm7KUMcA+7HvAIuAhwUsaOBO4/i8D95M6rJnDl13g6+bZydEmeD3B8HRugUqRLZCdSHnD8kHsouNIxgiSgSQsRw5XZocqQNwulbID2VtIs4HmwDbhSO4Fm8gNEADcbONBbETjQGwEO9CzIPldPyLsBqbxcYip67do10/Xr1023bt0y3bx5M/VzAgQpXr9+/Ux+McoDj7AlCkTwpLAKOADHOrA3LAzwlENTdQDzPIKjJODeQqLqAZfVEYR9aivVgt0wcMEuC9xDz4MdgRPhywBnbaQW8Fqw84ATU1Ae8JEjR7YRdRszZkwlYdkQAZiefO/ePdPdu3dN+gwLIU1jmAKwW7dupq5du7ZRBO5h+5U+BV1A0rDYGlbG+2kYLVwRqAFNMPbAy0CPoBsBLtCo04EjRkVUWeDybHJtZTfkw+0FLuXBbhdwAN+4ccN0+/Zt0/37902PHj0y8ffFixdTEPgUK32Alrp06VJTecCRgiNpXlyjwVIIklnA/aQqC3g96LVgdxpwwW4GcDKHPNgxC2Ek8Tf3e6hZayKIbbksNQqcLIj3MqIAH7OfFDgBMFrIs2fPTE+ePDHpGirOS/1aAqoHHFgo7imqx5YFLu9uFLik+8rAbhngPvUrAlvSypoHrokYeTWVYvOBQEbleL+AM9oagR3lN5vzxGwaa2HeQJmkmsAfPHhgevHihenDhw+mV69epXbDego+iUcxEYgTnHqwAS753fHOAJ4FutnAI3TVs13AAcC2WJwgoKLAEcGQSpA2+tw4blCTf0cAZWEXBc7uD+VhDYjOQDkk6mjAgYmVAPzNmzepmPQwGaJyRYHnwW4UeNaJgKLA68GWVK4ywP2imAddCPj79+9Nnz9/TqVeLpG5sP5MpQHOYZ14aEe+nwccYLIStuTYfQc0IPxZEaCzfMu9xALiDJkP9xSBXRQ475TygHvYBpxgSTrYLOC+VxFss8C3EvBYtjLAqQfWwoy3DfCHDx9WJLKRr1+/mn79+lX59OmT6fXr16Y7d+6YWFkEIgFL9hGHsIctCQZr1QgrYX07HsYhmPIMLR0z20XskscNZJ6t+7GdIrDrAfewI3B/TuWvBK7Pve9LtUDXgx5he3nYZFc0hgGPWcm3b99MP3/+rHz8+NEULQXwHPQEjKbQPoBo3UNwInCJAuozYMU9S/72h0slrc+zUcK9gIyNyPceOGWoBztLtWA3DFyw2wuchSa26poNnJ7Mu1nCZZSw/sIGit9EKQpc9lC018f0UOI91EPrPwlQv3//bhJoejlB8+XLlyZsBxti7VyVoEcBIg84YDTcyetZjOIZXMP32IZWCIHfkcA9tKJWE2H/VsDxX2yDzRGdewRePBrBKAE0J8MUa3w5JTYuCKKsNrLiWAQ4doGUCTEy887yiEPy48ePiuT/ffnyxfT27VvT06dPTQB//PixiRze74VyIBQwBFh/fkXCjvxh//8beIRdC7h6dEsC51QX4DlPzikBKZ4nxxYAzhKynsPzaTgshmt5vrcR3hPPwcdsqSxwiRmnOgqgYx15n5RgJUoFmfQQJAENYIIm6y+a7gOAaX9cfcwC7mEXBc5zygCnR3c0cEZmBM7+ASNMKg1cBWYzV8OSwgOcXo/icPbn0IHIvRQwa0Ob/3srkth+oxzAZThnicaiHsQdAfSBkHPxBGaJ75gr6H6eS7lRJnACowInFoJYvHr+/LmJHvI3Aye+MBp1Dc9lkY9RTv2oj+pvwNlsoEezrgJgggHBxa/ukX3g4UDj5bwUaZQwYti+Y+EMr2ZD2/cMGsn/JsnbQ/RNBSmGeN6eKMEz/hzFH33mu6YBF+DOAk4M6GjgZAVFgWeNlHhvPPWlDIlRBnDKwGjmWf4AbCIotACRPR7S5GW0qv8dI5/xcJ5Brg5coErRdgiwceWSoxlStB9/pt3L5/6kjoBFWAmd5R/wJgBXeX2skagPoMiGZAt8RvkZUUCM9/hz9XFpIv4SA+AKsol/KBfGwOSPmMXf+MSfyVFg7AF4bOFJ8UgGazPA9Uc0aDgai+/4G+uik/iKtgpwTa60Nt9pwD3sosDxdw9cmRIpK9cSoPROyhR/yo61+J+nSxqZWATlxw4YQfxNvf2PEriH5QqpFnAp0YtpeVoI32KYAkj/JwugMeLvNpngEDx1H0EQ6+J5NAC2Q/r57t27FCyzW7773YCzfNAywFkYk4CpZWIW1ZgTsAlCJoWlMDo0mvB34gm2QKPE7MevucS5ASOUevjAF0H79XC/gaLMjhkp6aQB52X+V8USBVQl4g88SRV5Od/zDJ4psPguvZJe62H/qcDVKTscuAoo3wWgLEICGhBZJEO6Ji6Y4eEIwNH/PXBGZjxz4y0LsATdOAfgGTSOTxW5hkZn5MYtQoD71PE/7eyQCLd/dVcAAAAASUVORK5CYII=",
[[43.60734015888762, 3.767730367838], [43.81775157192726, 4.3423021816889245]],
{"crossOrigin": false, "interactive": false, "opacity": 0.7, "zindex": 1}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
var layer_control_16772d432ce6954add8c65224c2554ae = {
base_layers : {
"openstreetmap" : tile_layer_8ebb7bdf73674872bc467e44a609f99a,
},
overlays : {
"January" : image_overlay_06ec3d8b1ee2472e0dee5d2551cc2776,
"February" : image_overlay_a89a77ffe85cc350b6d5e694476e9cb5,
"March" : image_overlay_9be1b876fc220c2b22bc75a44a970ca5,
"April" : image_overlay_d1707f042f199bc480911a410925d0ed,
"May" : image_overlay_047cf7a4389ddf495e7be954fc800a97,
"June" : image_overlay_35a70477cdb1aec56ce558f46abbb54b,
"July" : image_overlay_aa3f48d500ceccff387b5b78900cb83f,
"August" : image_overlay_6610ffc3843fc286b5f922d0ff999b64,
"September" : image_overlay_ba0f0d14f5a326f6d0321ceb4226f921,
},
};
L.control.layers(
layer_control_16772d432ce6954add8c65224c2554ae.base_layers,
layer_control_16772d432ce6954add8c65224c2554ae.overlays,
{"autoZIndex": true, "collapsed": true, "position": "topright"}
).addTo(map_c411caddeadbebf70f7cebfd52330cd5);
</script> </html>” width=”700” height=”500”style=”border:none !important;” “allowfullscreen” “webkitallowfullscreen” “mozallowfullscreen”></iframe>
Plot temporal series for the data
We calculate the mean temperature per month
data_interpolated
<xarray.DataArray 'LST_Day_1km' (time: 12, y: 48, x: 92)>
array([[[ 9.69, 9.69, 9.69, ..., 16.55, 16.71, 16.71],
[ 9.69, 9.69, 9.69, ..., 16.55, 16.71, 16.71],
[ 9.69, 9.69, 9.69, ..., 16.73, 16.69, 16.69],
...,
[ 15.51, 15.21, 15.21, ..., 11.71, 11.35, 11.35],
[ 15.47, 15.55, 15.55, ..., 11.71, 11.35, 11.35],
[ 15.47, 15.55, 15.55, ..., 13.47, 14.21, 14.21]],
[[ 12.93, 12.93, 12.93, ..., 17.77, 20.53, 20.53],
[ 12.93, 12.93, 12.93, ..., 20.53, 20.53, 20.53],
[ 12.93, 12.93, 12.93, ..., 20.53, 20.53, 20.53],
...,
[ 16.67, 17.83, 17.83, ..., 13.59, 14.23, 14.23],
[ 17.83, 17.83, 17.83, ..., 13.59, 14.23, 14.23],
[ 17.83, 17.83, 17.83, ..., 15.63, 15.97, 15.97]],
[[ 23.83, 23.27, 23.27, ..., 31.33, 31.65, 31.65],
[ 24.17, 23.93, 23.93, ..., 31.33, 31.65, 31.65],
[ 24.17, 23.93, 23.93, ..., 31.13, 30.55, 30.55],
...,
...
...,
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15],
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15],
[-273.15, -273.15, -273.15, ..., -273.15, -273.15, -273.15]],
[[ 8.65, 8.65, 8.65, ..., 9.33, 9.35, 9.35],
[ 8.63, 8.63, 8.63, ..., 9.33, 9.35, 9.35],
[ 8.63, 8.63, 8.63, ..., 9.55, 9.47, 9.47],
...,
[ 10.17, 10.23, 10.23, ..., 6.01, 6.09, 6.09],
[ 10.21, 10.17, 10.17, ..., 6.01, 6.09, 6.09],
[ 10.21, 10.17, 10.17, ..., 6.91, 7.67, 7.67]],
[[ 16.45, 16.07, 16.07, ..., 13.49, 13.49, 13.49],
[ 16.59, 16.17, 16.17, ..., 13.49, 13.49, 13.49],
[ 16.59, 16.17, 16.17, ..., 13.49, 13.49, 13.49],
...,
[ 16.65, 16.99, 16.99, ..., 12.91, 12.91, 12.91],
[ 16.63, 16.79, 16.79, ..., 12.91, 12.91, 12.91],
[ 16.63, 16.79, 16.79, ..., 12.91, 12.91, 12.91]]])
Coordinates:
* y (y) float64 6.303e+06 6.302e+06 ... 6.28e+06 6.279e+06
* x (x) float64 7.622e+05 7.628e+05 ... 8.072e+05 8.078e+05
* time (time) datetime64[ns] 2021-01-31 2021-02-28 ... 2021-12-31
spatial_ref int64 0
Attributes:
_FillValue: -273.15To plot temporal series of the temperature, we must calculate the averaage of temperatures per month in each raster. However, we need to deal with the nodata values first, since a high negative value can influence the mean value.
# Remove data where we have absolute kelvin -273.15. This will add np.nan where temperatures are -273.15
data_with_nan = [data_month.data for data_month in data_interpolated.where(data != -273.15)]
# Remove np.nan values
data_without_nan = [data_month[~np.isnan(data_month)] for data_month in data_with_nan]
# Calculate the mean with the present values
data_temp_mean = [np.nanmean(month_data) for month_data in data_without_nan ]
data_temp_mean
[11.963057065217424,
16.11516304347829,
26.789483695652205,
19.04933423913047,
29.070221920289885,
33.015901268115975,
39.149814311594234,
35.1516077898551,
19.180919384058,
-273.1500000000001,
9.299877717391333,
14.698691123188437]
The month of October has nodata values only. Lets replace it by the average of the temperatures of September and November
data_temp_mean[9] = np.mean( np.array([ data_temp_mean[8], data_temp_mean[10] ]), axis=0 )
data_temp_mean
[11.963057065217424,
16.11516304347829,
26.789483695652205,
19.04933423913047,
29.070221920289885,
33.015901268115975,
39.149814311594234,
35.1516077898551,
19.180919384058,
14.240398550724667,
9.299877717391333,
14.698691123188437]
We create a dataframe with the columns temperature and month
import pandas as pd
df_temp = pd.DataFrame({'month': list(months.keys()), 'temperature':data_temp_mean})
df_temp
| month | temperature | |
|---|---|---|
| 0 | January | 11.963057 |
| 1 | February | 16.115163 |
| 2 | March | 26.789484 |
| 3 | April | 19.049334 |
| 4 | May | 29.070222 |
| 5 | June | 33.015901 |
| 6 | July | 39.149814 |
| 7 | August | 35.151608 |
| 8 | September | 19.180919 |
| 9 | October | 14.240399 |
| 10 | november | 9.299878 |
| 11 | December | 14.698691 |
Plot data with seaborn
import seaborn as sns
sns.set_theme(rc={'figure.figsize':(8.7,7.27)})
sns.lineplot(x="month", y="temperature", data=df_temp)
plt.xticks(rotation=15)
plt.title(f'Montpellier verage temp in {year}')
plt.show()
Make a animation from the output raster
Create folders
print("Current Working Directory " , os.getcwd())
import matplotlib.pyplot as pyplot
in_dir = os.getcwd() + '/temp'
os.chdir(in_dir)
print("Input Directory " , in_dir)
Current Working Directory /content
Input Directory /content/temp
out_dir = os.getcwd() + '/animation/'
print("Animation out dir",out_dir)
if not os.path.exists(out_dir): os.makedirs(out_dir)
file_list = [file for file in os.listdir() if file.endswith('.tiff')]
print(os.listdir())
Animation out dir /content/temp/animation/
['modis_tile_10_.tiff', 'modis_tile_02_.tiff', 'modis_tile_03_.tiff', 'modis_tile_12_.tiff', 'modis_tile_05_.tiff', 'modis_tile_06_.tiff', 'modis_tile_11_.tiff', 'modis_tile_08_.tiff', 'animation', 'modis_tile_09_.tiff', 'modis_tile_04_.tiff', 'modis_tile_01_.tiff', 'modis_tile_07_.tiff']
for f in file_list:
raw_file = rasterio.open(f)
file = raw_file.read(1)
cmap='turbo'
plt.figure(figsize=(12,15))
plt.imshow(file,cmap=cmap,vmin=0,vmax=60)
parts = f.split("_")
file_name = parts[2]
print(f'Processing: {file_name}')
plt.title(f"Temperature {file_name}", fontsize=28,fontweight='bold')
cur_axes = plt.gca()
cur_axes.axes.get_xaxis().set_visible(False)
cur_axes.axes.get_yaxis().set_visible(False)
cb = plt.colorbar(fraction=0.046, pad=0.04,shrink=0.3)
cb.set_label('Temperature time series', size=24)
cb.ax.tick_params(labelsize=18)
plt.tight_layout()
print(out_dir)
plt.savefig('{}{}_temperature.png'.format(out_dir,file_name,dpi=150))
plt.close()
Processing: 10
/content/temp/animation/
Processing: 02
/content/temp/animation/
Processing: 03
/content/temp/animation/
Processing: 12
/content/temp/animation/
Processing: 05
/content/temp/animation/
Processing: 06
/content/temp/animation/
Processing: 11
/content/temp/animation/
Processing: 08
/content/temp/animation/
Processing: 09
/content/temp/animation/
Processing: 04
/content/temp/animation/
Processing: 01
/content/temp/animation/
Processing: 07
/content/temp/animation/
Create animation
import imageio.v2 as imageio
def make_gif(input_folder, save_file_path):
episode_frames = []
time_per_step = 4
for root, _, files in os.walk(input_folder):
file_paths = [os.path.join(root,file) for file in files]
file_paths = sorted (file_paths,key=lambda x:os.path.getmtime(x))
episode_frames = [imageio.imread(file_path) for file_path in file_paths if file_path.endswith('.png')]
episode_frames = np.array(episode_frames)
imageio.mimsave(save_file_path,episode_frames,duration=6)
make_gif(out_dir,os.path.join(out_dir,"Animation_Temp.gif"))
Visualize the animation
from IPython.display import Image, display
with open(out_dir + "Animation_Temp.gif", 'rb') as f:
display(Image(data=f.read(),format='png',width = 600, height = 600))
