diff --git a/scripts/ add_heat_data/ add_heat_data.py b/scripts/ add_heat_data/ add_heat_data.py new file mode 100644 index 0000000..82ab550 --- /dev/null +++ b/scripts/ add_heat_data/ add_heat_data.py @@ -0,0 +1,200 @@ +import ssl +import urllib.request +import json +import rasterio +import geopandas as gpd +from rasterio.mask import mask +import numpy as np +from affine import Affine +from math import cos, radians +import psycopg2 +from psycopg2.extras import RealDictCursor +from shapely.geometry import shape +from dotenv import dotenv_values +from math import cos, radians + +# Load environment variables from .env file +env_vars = dotenv_values(".env") +date = env_vars["DATE"] +# Read min and max temperature from environment variables +min_temp_k = float(env_vars["MIN_TEMP_K"]) +max_temp_k = float(env_vars["MAX_TEMP_K"]) + +def fetch_wfs_data(typename): + + # Base URL for the WFS server + wfs_url = 'https://kartta.hsy.fi/geoserver/wfs' + + # Parameters for the WFS request + params = { + 'service': 'WFS', + 'request': 'GetFeature', + 'typename': typename, + 'version': '2.0.0', + 'outputFormat': 'application/json', + 'srsName': 'EPSG:4326' + } + + # Construct the query string + query_string = urllib.parse.urlencode(params) + + # Complete URL for the request + request_url = wfs_url + '?' + query_string + + print(request_url) + + # Create an SSL context that does not verify certificates + ssl_context = ssl.create_default_context() + ssl_context.check_hostname = False + ssl_context.verify_mode = ssl.CERT_NONE + + # Request to WFS and process response + try: + response = urllib.request.urlopen(request_url, context=ssl_context) + data = json.loads(response.read()) + return data + except Exception as e: + print(f"An error occurred: {e}") + +def calculate_cell_area(transform, latitude): + """ + Calculate the area of a cell in square meters for a raster in geographic coordinates, + assuming a uniform cell size and using an approximate method suitable for the latitude of Helsinki. + + Parameters: + - transform: The affine transform of the raster. + - latitude: The latitude at which to calculate the area (default is Helsinki's approximate latitude). + + Returns: + - The area of a cell in square meters. + """ + # Calculate dimensions of a cell in degrees + delta_lat = abs(transform.e) # Cell height in degrees (transform.e is negative) + delta_lon = transform.a # Cell width in degrees + + # Convert latitude to radians for the cosine calculation + lat_rad = radians(latitude) + + # Calculate the area of a cell using the approximate formula + cell_area_km2 = (delta_lon * cos(lat_rad)) * delta_lat * ((40075.017/360) ** 2) # Earth circumference is about 40075 km + + # Convert square kilometers to square meters + cell_area_m2 = cell_area_km2 * (1000 ** 2) + + return cell_area_m2 + +def calculate_heat_exp(geojson_data, raster_path, output_column_name): + # Open the raster and vector files + gdf = gpd.GeoDataFrame.from_features(geojson_data["features"]) + src = rasterio.open(raster_path) + + # Initialize a list to store weighted average values + weighted_avg_values = [] + + # Loop through each building polygon + for idx, geom in gdf.iterrows(): + # Extract the centroid of the geometry to use its latitude + centroid = geom.geometry.centroid + latitude = centroid.y + geometry = geom.geometry.__geo_interface__ + out_image, out_transform = mask(src, [geometry], crop=True, all_touched=True, filled=True) + out_image = out_image[0] # Assuming a single band + + # Create an Affine transform for the masked region + masked_transform = Affine(out_transform[0], out_transform[1], out_transform[2],out_transform[3], out_transform[4], out_transform[5]) + + # Initialize total value and total area + total_value = 0 + total_area = 0 + + # Calculate total value and total area for the masked region + for row in range(out_image.shape[0]): + for col in range(out_image.shape[1]): + pixel_value = out_image[row, col] + if pixel_value > 0: # Valid pixel + cell_area = calculate_cell_area(masked_transform, latitude) + total_value += pixel_value * cell_area + total_area += cell_area + + # Calculate the weighted average value if there's any valid area + if total_area > 0: + avg_value = total_value / total_area + else: + avg_value = np.nan # No valid data within this polygon + + weighted_avg_values.append(avg_value) + + # Add weighted average heat values to the GeoDataFrame + gdf[output_column_name] = weighted_avg_values + + return gdf + +def insert_to_db(buildings_gdf): + """ + Compare fetched WFS data with postal code GeoJSON file, join posno, and save to PostgreSQL table. + + Parameters: + - wfs_json (dict): Fetched data from the WFS server. + - postal_code_geojson_path (str): Path to the postal code GeoJSON file. + """ + # Read postal code data + + # Read environment variables (ensure they are set properly) + DB_HOST = env_vars["DB_HOST"] + DB_PORT = env_vars["DB_PORT"] + DB_NAME = env_vars["DB_NAME"] + DB_USER = env_vars["DB_USER"] + DB_PASSWORD = env_vars["DB_PASSWORD"] + HEAT_TIMESERIES_TABLE = env_vars["HEAT_TIMESERIES_TABLE"] + + if None in [DB_HOST, DB_PORT, DB_NAME, DB_USER, DB_PASSWORD, HEAT_TIMESERIES_TABLE]: + raise ValueError("One or more required environment variables are missing.") + + + insert_query = f""" + INSERT INTO {HEAT_TIMESERIES_TABLE} (avgheatexposure, date, vtj_prt, avg_temp_c) + VALUES (%s, %s, %s, %s) + """ + + # Prepare data for insertion + data_to_insert = [ + (row.avgheatexposure, date, row.vtj_prt, calculate_temp_in_c(row.avgheatexposure)) + for _, row in buildings_gdf.iterrows() if row.avgheatexposure >= 0 + ] + + # Insert data into database + try: + with psycopg2.connect( + host=DB_HOST, port=DB_PORT, database=DB_NAME, + user=DB_USER, password=DB_PASSWORD + ) as conn: + with conn.cursor() as cur: + cur.executemany(insert_query, data_to_insert) + conn.commit() + except Exception as e: + print(f"Error inserting data into database: {e}") + raise + +def calculate_temp_in_c(heatexposure): + """ + Calculates temperature in Celsius from a normalized heatexposure index and reference temperatures in Kelvin. + + Args: + heatexposure (float): The normalized heatexposure index. + max_temp_k (float): The maximum reference temperature in Kelvin. + min_temp_k (float): The minimum reference temperature in Kelvin. + + Returns: + float: The calculated temperature in Celsius. + """ + temp_k = heatexposure * (max_temp_k - min_temp_k) + min_temp_k # Denormalize to Kelvin + temp_c = temp_k - 273.15 # Convert from Kelvin to Celsius + return temp_c + +typename = 'asuminen_ja_maankaytto:pks_rakennukset_paivittyva' +raster_path = env_vars["RASTER"] +wfs_json = fetch_wfs_data(typename) +gdf = calculate_heat_exp(wfs_json, raster_path, 'avgheatexposure' ) +print(f"Number of filtered features: {len(gdf)}") + +insert_to_db( gdf ) \ No newline at end of file diff --git a/scripts/ add_heat_data/requirements.txt b/scripts/ add_heat_data/requirements.txt new file mode 100644 index 0000000..0a0b8fb --- /dev/null +++ b/scripts/ add_heat_data/requirements.txt @@ -0,0 +1,8 @@ +geopandas +rasterio +psycopg2 +shapely +python-dotenv +numpy +affine +webcolors \ No newline at end of file diff --git a/scripts/heat_exposure/calculate_and_insert.py b/scripts/heat_exposure/calculate_and_insert.py index 0ea3ece..82ab550 100644 --- a/scripts/heat_exposure/calculate_and_insert.py +++ b/scripts/heat_exposure/calculate_and_insert.py @@ -192,7 +192,7 @@ def calculate_temp_in_c(heatexposure): return temp_c typename = 'asuminen_ja_maankaytto:pks_rakennukset_paivittyva' -raster_path = 'normalised.tiff' +raster_path = env_vars["RASTER"] wfs_json = fetch_wfs_data(typename) gdf = calculate_heat_exp(wfs_json, raster_path, 'avgheatexposure' ) print(f"Number of filtered features: {len(gdf)}") diff --git a/scripts/hsylandcover.py b/scripts/hsylandcover/hsylandcover.py similarity index 100% rename from scripts/hsylandcover.py rename to scripts/hsylandcover/hsylandcover.py diff --git a/scripts/hsylandcover/requirements.txt b/scripts/hsylandcover/requirements.txt new file mode 100644 index 0000000..f9b0b84 --- /dev/null +++ b/scripts/hsylandcover/requirements.txt @@ -0,0 +1,3 @@ +google-cloud-sql +geopandas +pandas \ No newline at end of file diff --git a/scripts/hsytrees.py b/scripts/hsytrees.py deleted file mode 100644 index caf86fc..0000000 --- a/scripts/hsytrees.py +++ /dev/null @@ -1,122 +0,0 @@ -import ssl -import urllib.request -import json -import geopandas as gpd -import psycopg2 -from psycopg2.extras import RealDictCursor -from shapely.geometry import shape -from dotenv import dotenv_values -import os - -# Load environment variables from .env file -env_vars = dotenv_values(".env") - -def fetch_wfs_data(typename, city): - - # Base URL for the WFS server - wfs_url = 'https://kartta.hsy.fi/geoserver/wfs' - - # Parameters for the WFS request - params = { - 'service': 'WFS', - 'request': 'GetFeature', - 'typename': typename, - 'version': '2.0.0', - 'outputFormat': 'application/json', - 'CQL_FILTER': 'kunta IN ' + city, - 'srsName': 'EPSG:4326' - } - - # Construct the query string - query_string = urllib.parse.urlencode(params) - - # Complete URL for the request - request_url = wfs_url + '?' + query_string - - print(request_url) - - # Create an SSL context that does not verify certificates - ssl_context = ssl.create_default_context() - ssl_context.check_hostname = False - ssl_context.verify_mode = ssl.CERT_NONE - - # Request to WFS and process response - try: - response = urllib.request.urlopen(request_url, context=ssl_context) - data = json.loads(response.read()) - return data - except Exception as e: - print(f"An error occurred: {e}") - -def spatial_relation_with_postal_code_areas_and_save_to_db(wfs_json, postal_code_geojson_path): - """ - Compare fetched WFS data with postal code GeoJSON file, join posno, and save to PostgreSQL table. - - Parameters: - - wfs_json (dict): Fetched data from the WFS server. - - postal_code_geojson_path (str): Path to the postal code GeoJSON file. - """ - # Read postal code data - - conn = psycopg2.connect( - host=env_vars["DB_HOST"], - port=env_vars["DB_PORT"], - database=env_vars["DB_NAME"], - user=env_vars["DB_USER"], - password=env_vars["DB_PASSWORD"] - ) - cur = conn.cursor() - - # Get table name from environment variables - table_name = env_vars["TABLE_NAME"] - - postal_code_gdf = gpd.read_file(postal_code_geojson_path) - wfs_gdf = gpd.GeoDataFrame.from_features(wfs_json["features"]) - print(f"Number of filtered features: {len(wfs_gdf)}") - - wfs_gdf.crs = "EPSG:4326" - postal_code_gdf.crs = "EPSG:4326" - - # Iterate over each feature in postal_code_gdf - for _, postalarea in postal_code_gdf.iterrows(): - try: - # Clip the hsylandcover features by the district geometry - clipped = wfs_gdf.clip(postalarea.geometry) - - # Skip if there are no clipped geometries - if clipped.empty: - continue - - # Calculate area for each clipped feature and sum up - for _, clipped_feature in clipped.iterrows(): - # Get the geometry in WKT format - if clipped_feature['geometry'].geom_type == 'Polygon': - geometry_wkt = clipped_feature['geometry'].wkt - elif clipped_feature['geometry'].geom_type == 'MultiPolygon': - geometry_wkt = clipped_feature['geometry'].buffer(0).wkt # Handle invalid geometries - else: - print(f"Unsupported geometry type: {clipped_feature['geometry'].geom_type}") - continue - - # Insert to database - cur.execute(""" - INSERT INTO {table_name} (postinumero, kohde_id, kunta, paaluokka, alaluokka, ryhma, koodi, kuvaus, geom) - VALUES (%s, %s, %s, %s, %s, %s, %s, %s, ST_GeomFromText(%s)) - """, (postalarea['posno'], clipped_feature['kohde_id'], clipped_feature['kunta'], clipped_feature['paaluokka'], clipped_feature['alaluokka'], clipped_feature['ryhma'], clipped_feature['koodi'], clipped_feature['kuvaus'], geometry_wkt)) - - # Commit after processing each postal area - conn.commit() - print(f"Data saved to PostgreSQL table tree_f for postal area {postalarea['posno']} successfully!") - - except Exception as e: - print(f"An error occurred while processing postalarea {postalarea['posno']}: {e}") - - # Close the cursor and the connection - cur.close() - conn.close() - -typename = env_vars["TYPENAME"] -city = env_vars["CITY"] -postal_code_geojson_path = 'hsy_po.json' -wfs_json = fetch_wfs_data(typename, city) -spatial_relation_with_postal_code_areas_and_save_to_db(wfs_json, postal_code_geojson_path) \ No newline at end of file diff --git a/scripts/hsytrees/hsytrees.py b/scripts/hsytrees/hsytrees.py new file mode 100644 index 0000000..abd323e --- /dev/null +++ b/scripts/hsytrees/hsytrees.py @@ -0,0 +1,111 @@ +import os +import json +import urllib.request +import urllib.parse +import ssl +import geopandas as gpd +from shapely.geometry import shape +from google.cloud import sql_v1beta4 + +def fetch_wfs_data(typename, city): + + # Base URL for the WFS server + wfs_url = 'https://kartta.hsy.fi/geoserver/wfs' + + # Parameters for the WFS request + params = { + 'service': 'WFS', + 'request': 'GetFeature', + 'typename': typename, + 'version': '2.0.0', + 'outputFormat': 'application/json', + 'CQL_FILTER': 'kunta IN ' + city, + 'srsName': 'EPSG:4326' + } + + # Construct the query string + query_string = urllib.parse.urlencode(params) + + # Complete URL for the request + request_url = wfs_url + '?' + query_string + + # Create an SSL context that does not verify certificates + ssl_context = ssl.create_default_context() + ssl_context.check_hostname = False + ssl_context.verify_mode = ssl.CERT_NONE + + # Request to WFS and process response + try: + response = urllib.request.urlopen(request_url, context=ssl_context) + data = json.loads(response.read()) + return data + except Exception as e: + print(f"An error occurred: {e}") + +def spatial_relation_with_postal_code_areas_and_save_to_db(wfs_json, postal_code_geojson_path): + """ + Compare fetched WFS data with postal code GeoJSON file, join posno, and save to PostgreSQL table. + """ + + # Google Cloud SQL connection details + instance_name = os.environ.get("INSTANCE_NAME") + db_name = os.environ.get("DB_NAME") + db_user = os.environ.get("DB_USER") + db_password = os.environ.get("DB_PASSWORD") + + db_config = { + "user": db_user, + "password": db_password, + "database": db_name, + } + pool = sql_v1beta4.ConnectionPool(instance_connection_name=instance_name, db_config=db_config) + table_name = os.environ.get("TABLE_NAME") # Get table name from environment variables + + postal_code_gdf = gpd.read_file(postal_code_geojson_path) + wfs_gdf = gpd.GeoDataFrame.from_features(wfs_json["features"]) + print(f"Number of filtered features: {len(wfs_gdf)}") + + wfs_gdf.crs = "EPSG:4326" + postal_code_gdf.crs = "EPSG:4326" + + with pool.connect() as conn: + with conn.cursor() as cursor: + # Iterate over each feature in postal_code_gdf + for _, postalarea in postal_code_gdf.iterrows(): + try: + # Clip the hsylandcover features by the district geometry + clipped = wfs_gdf.clip(postalarea.geometry) + + # Skip if there are no clipped geometries + if clipped.empty: + continue + + # Calculate area for each clipped feature and sum up + for _, clipped_feature in clipped.iterrows(): + # Get the geometry in WKT format + if clipped_feature['geometry'].geom_type == 'Polygon': + geometry_wkt = clipped_feature['geometry'].wkt + elif clipped_feature['geometry'].geom_type == 'MultiPolygon': + geometry_wkt = clipped_feature['geometry'].buffer(0).wkt # Handle invalid geometries + else: + print(f"Unsupported geometry type: {clipped_feature['geometry'].geom_type}") + continue + + # Insert to database + cursor.execute(f""" + INSERT INTO {table_name} (postinumero, kohde_id, kunta, paaluokka, alaluokka, ryhma, koodi, kuvaus, geom) + VALUES (%s, %s, %s, %s, %s, %s, %s, %s, ST_GeomFromText(%s)) + """, (postalarea['posno'], clipped_feature['kohde_id'], clipped_feature['kunta'], clipped_feature['paaluokka'], clipped_feature['alaluokka'], clipped_feature['ryhma'], clipped_feature['koodi'], clipped_feature['kuvaus'], geometry_wkt)) + + # Commit after processing each postal area + conn.commit() + print(f"Data saved to PostgreSQL table {table_name} for postal area {postalarea['posno']} successfully!") + + except Exception as e: + print(f"An error occurred while processing postalarea {postalarea['posno']}: {e}") + +typename = os.environ.get["TYPENAME"] +city = os.environ.get["CITY"] +postal_code_geojson_path = os.environ.get["hsy_po"] +wfs_json = fetch_wfs_data(typename, city) +spatial_relation_with_postal_code_areas_and_save_to_db(wfs_json, postal_code_geojson_path) \ No newline at end of file diff --git a/scripts/hsytrees/requirements.txt b/scripts/hsytrees/requirements.txt new file mode 100644 index 0000000..26132c6 --- /dev/null +++ b/scripts/hsytrees/requirements.txt @@ -0,0 +1,4 @@ +google-cloud-sql +geopandas +shapely +psycopg2-binary \ No newline at end of file