<jupyter_start><jupyter_text>Une introduction pragmatique à SQL pour la géomatique--- Cours 5 : SQL raster---  Installation de la librairie *rasterio* pour traiter des images raster en python<jupyter_code>#!pip install rasterio --user
# Librairies python
import geopandas as gpd
import getpass
from matplotlib import pyplot
import psycopg2
from rasterio.io import MemoryFile
from rasterio.plot import show<jupyter_output><empty_output><jupyter_text>Connexion à PostGIS<jupyter_code>ip = '192.168.3.11'
uid = input('Entrez votre indentifiant :')
mdp = getpass.getpass('Entrez votre mot de passe :')
bd = 'climate'
%load_ext sql
%sql postgresql://$uid:$mdp@$ip/$bd<jupyter_output><empty_output><jupyter_text>*** World ClimateLes images raster (.tif, .asc etc.) sont préparées et chargées dans PostGIS grâce à des outils comme [raster2pgsql](https://postgis.net/docs/using_raster_dataman.htmlRT_Raster_Loader) et/ou [gdal](https://gdal.org). Dans ce cours nous nous concentrerons sur les requêtes SQL.Les jeux de données que nous allons utilisés proviennent de la base de données "WorldClim". Par exemple, la moyenne des températures minimales sur la période 1970 - 2000. Ici chaque pixel représente environ 1 $km^2$. La valeur du pixel donne cette température car il s'agit d'un raster mono-bande. Dans un raster multi-bande, à chaque pixel sont associées plusieurs valeurs représentant différentes mesures (températures, précipitations, couleurs etc.)  Dans PostGIS, chaque enregistrement de une table raster est une dalle de pixels.<jupyter_code>conn = psycopg2.connect("postgresql://" + uid + ":" +  mdp + "@" + ip + "/" + bd)
cur = conn.cursor()
# ensure that the GTiff driver is available, 
# see https://postgis.net/docs/postgis_gdal_enabled_drivers.html
cur.execute('''
  SET postgis.gdal_enabled_drivers TO 'GTiff';
  WITH temp as (
    SELECT ST_Union(rast,1) As file_rast
    FROM wc21_30s_bio_6_france
  )
  SELECT ST_AsGDALRaster(file_rast, 'GTiff') FROM temp;
 ''')

(fig, ax) = pyplot.subplots(1, 1, figsize = (10, 10))
for row in cur:
    rast = row[0].tobytes()
    with MemoryFile(rast).open() as dataset:
        show(dataset.read(1), cmap = 'Greens', ax = ax, alpha = 0.7, title = 'Emprise du raster')
%%sql
SELECT st_height(rast) as hauteur, st_width(rast) as largeur
FROM wc21_30s_bio_6_france
WHERE rid = 1; -- identifiant de la dalle<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
1 rows affected.<jupyter_text>Température au centroïde de Grenoble :<jupyter_code>%%sql
SELECT nom_com_m, ST_value(rast, st_centroid(geom)) as temp_min
FROM commune, wc21_30s_bio_6_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, st_centroid(geom));<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
1 rows affected.<jupyter_text>On peut avoir la valeur de n'importe quel pixel de la dalle en faisant varier les 2 derniers paramètres entre 86 et 81 :<jupyter_code>%%sql
SELECT nom_com_m, ST_value(rast, 80, 85) as temp_min
FROM commune, wc21_30s_bio_6_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, st_centroid(geom));<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
1 rows affected.<jupyter_text>On peut utiliser des objets vectoriels pour géotraiter les raster. Par exemple en utilisant la fonction st_intersection(). Les objets obtenus sont de type 'geomval' (ici polygones représentant une même température) dont les valeurs sont accessibles par la notation pointée .val :<jupyter_code>%%sql

-- Températures minimales à Grenoble

SELECT nom_com_m, st_astext((st_intersection(rast, geom)).geom) as geom, (st_intersection(rast, geom)).val
FROM commune, wc21_30s_bio_6_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, geom)
LIMIT 2;<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
2 rows affected.<jupyter_text>La température minimale moyenne à Grenoble est obtenue par :<jupyter_code>%%sql
-- Température minimale moyenne à Grenoble

DROP VIEW IF EXISTS decoupage_grenoble;

CREATE VIEW decoupage_grenoble AS
SELECT nom_com_m, (st_intersection(rast, geom)).val
FROM commune, wc21_30s_bio_6_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, geom);

SELECT nom_com_m, avg(val) AS moyenne_temp_min
FROM decoupage_grenoble
GROUP BY nom_com_m;<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
Done.
Done.
1 rows affected.<jupyter_text>De façon équivalente :<jupyter_code>%%sql

-- Température minimale moyenne à Grenoble

WITH decoupage_grenoble AS (
SELECT nom_com_m, (st_intersection(rast, geom)).val
FROM commune, wc21_30s_bio_6_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, geom)
)
SELECT nom_com_m, avg(val) AS moyenne_temp_min
FROM decoupage_grenoble
GROUP BY nom_com_m;
%%sql

-- Température maximale moyenne à Grenoble

WITH decoupage_grenoble AS (
SELECT nom_com_m, (st_intersection(rast, geom)).val
FROM commune, wc21_30s_bio_5_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, geom)
)
SELECT nom_com_m, avg(val) AS moyenne_temp_max
FROM decoupage_grenoble
GROUP BY nom_com_m;
%%sql

-- Communes les plus froides
WITH decoupage_metro AS (
SELECT nom, (st_intersection(rast, geom)).val
FROM limites_communales_metro_epgs4326, wc21_30s_bio_6_france
WHERE st_intersects(rast, geom)
)
SELECT nom, avg(val) AS moyenne_temp_min
FROM decoupage_metro
GROUP BY nom
ORDER BY moyenne_temp_min ASC 
LIMIT 5;
%%sql

-- Communes les plus chaudes
WITH valeur_des_differents_raster AS (
SELECT nom, (st_intersection(rast, geom)).val
FROM limites_communales_metro_epgs4326, wc21_30s_bio_5_france
WHERE st_intersects(rast, geom)
)
SELECT nom, avg(val) AS moyenne_temp_max
FROM valeur_des_differents_raster
GROUP BY nom
ORDER BY moyenne_temp_max DESC 
LIMIT 5;
%%sql

-- Communes les plus élevées
WITH decoupage_metro AS (
SELECT nom, (st_intersection(rast, geom)).val
FROM limites_communales_metro_epgs4326, wc21_30s_elev_france
WHERE st_intersects(rast, geom)
)
SELECT nom, avg(val) AS moyenne_altitude
FROM decoupage_metro
GROUP BY nom
ORDER BY moyenne_altitude DESC 
LIMIT 5;
%%sql

-- Communes les plus plates
WITH valeur_des_differents_raster AS (
SELECT nom, (st_intersection(rast, geom)).val
FROM limites_communales_metro_epgs4326, wc21_30s_elev_france
WHERE st_intersects(rast, geom)
)
SELECT nom, stddev(val) AS stddev_altitude
FROM valeur_des_differents_raster
GROUP BY nom
ORDER BY stddev_altitude ASC 
LIMIT 5;<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
5 rows affected.<jupyter_text>Des modèles climatiques (par exemple CNRM-CM6-1 de Météo France) permettent d'estimer l'évolution du climat selon différents scénarios socio-économiques (Shared Socio-economic Pathways, SSPs: 126, 245, 370 et 585).<jupyter_code>%%html
<iframe scrolling="no" frameborder="0" marginheight="0px" marginwidth="0px" style="display: initial; margin: 0 auto;" src="https://cbhighcharts2019.s3.eu-west-2.amazonaws.com/CMIP6/emissions+cmip6.html" width="770px" height="500px"></iframe><span style="display:block; height:22px; max-width:800px;"><a href="https://www.carbonbrief.org"><img src="https://s3.eu-west-2.amazonaws.com/cbhighcharts2019/cb-logo-highcharts.svg" style="width: 22px; height: 22px; margin-top: 2px; margin-bottom: 2px; float:right; background-repeat: no-repeat; background-size: contain;"/></a></span><jupyter_output><empty_output><jupyter_text>Sur la période 2020 - 2040, quelles sont les prévisions pour Grenoble ?<jupyter_code>%%sql

-- Température maximale moyenne à Grenoble
-- wc21_25m_bioc_cnrm_cm61_ssp126_2021_2040_france est un raster multi-bande
-- Chacune des 19 bandes correspond à l'une des 19 variables bioclimatiques
-- bio5 = tmax
-- bio6 = tmin

WITH decoupage_grenoble AS (
SELECT nom_com_m, (st_intersection(geom, rast, 5)).val
FROM commune, wc21_25m_bioc_cnrm_cm61_ssp126_2021_2040_france
WHERE nom_com_m = 'GRENOBLE' AND st_intersects(rast, geom)
)
SELECT nom_com_m, avg(val) AS temp_max_ssp126
FROM decoupage_grenoble
GROUP BY nom_com_m;<jupyter_output>postgresql://fiorinhu:***@192.168.3.11/bd_fiorinhu
 * postgresql://fiorinhu:***@192.168.3.11/climate
1 rows affected.<jupyter_text>**Exercice :** déterminez t_max pour les autres scénarios. Faites de même pour t_min.  *** Raz de maréeSur la Côte d'Azur, les simulations prévoient des vagues jusqu'à 4 m en cas de séisme (la montée des eaux due au changement climatique est estimée à 1 m).  Voyons quel pourrait être l'impact (en l'absence de tout autre considération) d'une telle vague sur Saint Laurent-du-Var et Cagnes-sur-mer.<jupyter_code>conn = psycopg2.connect("postgresql://" + uid + ":" +  mdp + "@" + ip + "/" + bd)
cur = conn.cursor()

bati = gpd.read_postgis("SELECT * FROM nice.bati_cagnes_st_laurent;", con = conn)

# ensure that the GTiff driver is available, 
# see https://postgis.net/docs/postgis_gdal_enabled_drivers.html
cur.execute('''
  SET postgis.gdal_enabled_drivers TO 'GTiff';
  WITH temp as (
    SELECT ST_Union(rast,1) As file_rast
    FROM nice.o_4_cagnes_st_laurent
  )
  SELECT ST_AsGDALRaster(file_rast, 'GTiff') FROM temp;
 ''')
(fig, (axl, axr)) = pyplot.subplots(1, 2, figsize = (21, 15))
bati.plot(ax = axr)
pyplot.title('bâtiments')
axr.set_axis_off()
for row in cur:
    rast = row[0].tobytes()
    with MemoryFile(rast).open() as dataset:
        show(dataset.read(1), cmap = 'Greens', ax = axl, alpha = 0.7, title = 'mnt')<jupyter_output><empty_output><jupyter_text>Combien de bâtiments sont situés entre 0 et 4 m ?<jupyter_code>%%sql

SELECT count(*) as nb_batiments
FROM nice.bati_cagnes_st_laurent, nice.cagnes_st_laurent -- RGEALTI 1m
WHERE ST_Intersects(geom, rast) AND ST_Value(rast, ST_Centroid(geom)) < 5;
sql = '''
SELECT geom 
FROM nice.bati_cagnes_st_laurent, nice.cagnes_st_laurent 
WHERE ST_Intersects(geom, rast) AND ST_Value(rast, ST_Centroid(geom)) < 5;
'''
bati_0_5 = gpd.read_postgis(sql, con = conn)
bati_0_5.plot().set_axis_off()<jupyter_output><empty_output><jupyter_text>Combien de personnes vivent dans cette frange littorale (hors tourisme) ?  Pour cela, nous allons utiliser les données carroyées 200 m de l'INSEE.<jupyter_code>%%sql
WITH bande_0_5 as (
SELECT geom 
FROM nice.bati_cagnes_st_laurent, nice.cagnes_st_laurent 
WHERE ST_Intersects(geom, rast) AND ST_Value(rast, ST_Centroid(geom)) < 5
)
SELECT  sum(ind) as personnes, 
        sum(ind_65_79 + ind_80p) as plus_de_65, 
        sum(log_ap90) as logements_apres_90,
        ceil(sum(ind_snv)/sum(men)) as niveau_de_vie_par_menage
FROM nice.carreaux
WHERE idinspire IN (SELECT DISTINCT idinspire
                    FROM bande_0_5 as B, nice.carreaux as C
                    WHERE ST_Intersects(B.geom, C.geom));<jupyter_output>* postgresql://fiorinhu:***@192.168.3.11/climate
1 rows affected.
{'fraction_non_alphanumeric': 0.10742596810933941, 'fraction_numerical': 0.03744874715261959, 'mean_word_length': 5.792079207920792, 'pattern_counts': {'":': 2, '<': 55, '<?xml version=': 0, '>': 52, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 16, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '86660', 'license': 'CC0-1.0', 'n_tokens_mistral': 4558, 'n_tokens_neox': 4244, 'n_words': 1025, 'path': 'postgis/sql_5.ipynb', 'repo_name': 'fiorinoh/geomas-cours', 'size': '11007'}