turns-00029.parquet:35049
1c17b9d3973d4c008d06abe8
turn 1/1o1-preview-2024-09-12EnglishPoland6940 words
degenerate_repetitionAbsentFinal dense release
USER
import geopandas as gpd
import matplotlib.pyplot as plt
from geopy.geocoders import Nominatim
import requests
from bs4 import BeautifulSoup
import re
from tqdm import tqdm
from shapely.geometry import Point
import json
import pandas as pd
import numpy as np
import rasterio
import warnings
from scipy.spatial import cKDTree
from rasterio.transform import from_origin
from rasterio.features import rasterize
from rasterio.warp import reproject, Resampling
import os
import rasterio.mask
import matplotlib.colors as mcolors
import matplotlib.patheffects as patheffects
from matplotlib.font_manager import FontProperties
from pykrige.uk import UniversalKriging
from sklearn.linear_model import HuberRegressor
from scipy.stats import zscore
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from scipy.ndimage import gaussian_filter
from datetime import datetime, time
import pytz
import math
import glob
from rasterio.merge import merge
os.environ['PROJ_LIB'] = r'C:\Users\Lenovo\Documents\conda\Library\share\proj'
warnings.filterwarnings("ignore")
# Keywords to remove
keywords_to_remove = [
r"(centrum)", "rondo", "most", "Obwodnica", "obwodnica I", "obwodnica II", "obwodnica",
"- Zawiszyn", "OUD Węzeł", "PPO", "Węzeł", "N (wąwóz)", "S (most)", "stacja", "z pomiarów", "min", "max", "średnia"
]
def clean_station_name(station_name):
pattern = re.compile("|".join(keywords_to_remove), re.IGNORECASE)
cleaned_name = pattern.sub("", station_name).strip()
return cleaned_name
def get_masovian_boundary():
url = "https://raw.githubusercontent.com/ppatrzyk/polska-geojson/master/wojewodztwa/wojewodztwa-medium.geojson"
# Fetch the content using requests
response = requests.get(url)
if response.status_code == 200:
data = response.json()
# Set the CRS to EPSG:4326 when creating the GeoDataFrame
gdf = gpd.GeoDataFrame.from_features(data["features"], crs="EPSG:4326")
masovian = gdf[gdf['nazwa'] == 'mazowieckie']
return masovian
else:
print(f"Error fetching GeoJSON data: {response.status_code}")
return None
# Global variable for Masovian boundary
masovian_boundary = get_masovian_boundary()
def is_in_masovian(lat, lon):
point = gpd.GeoSeries([Point(lon, lat)], crs="EPSG:4326")
masovian_geom = masovian_boundary.to_crs(epsg=2180).reset_index(drop=True)
point_projected = point.to_crs(epsg=2180)
return masovian_geom.contains(point_projected.buffer(100)).values[0]
def fetch_traxelektronik_data():
url = 'https://www.traxelektronik.pl/pogoda/zbiorcza.php?RejID=10'
response = requests.get(url)
soup = BeautifulSoup(response.text, 'html.parser')
stations = []
temperatures = []
table_rows = soup.find_all('tr')[2:]
for row in table_rows:
columns = row.find_all('td')
if len(columns) > 1:
station = columns[0].get_text(strip=True)
temperature = columns[1].get_text(strip=True)
if temperature != "-":
stations.append(station) # Keep original station names
temperatures.append(temperature)
return stations, temperatures
def fetch_netatmo_data():
url = "https://api.netatmo.com/api/getpublicdata"
params = {
"lat_ne": 53,
"lon_ne": 22,
"lat_sw": 51,
"lon_sw": 19,
"required_data": "temperature",
"filter": "true",
"access_token": "YOUR_ACCESS_TOKEN"
}
response = requests.get(url, params=params)
if response.status_code != 200:
print(f"Error fetching Netatmo data: {response.status_code}")
return [], [], []
data = json.loads(response.text)
stations = []
temperatures = []
coordinates = []
# List of cities to exclude
excluded_cities = ["Płońsk", "Łochów"]
for station in data.get('body', []):
place = station.get('place', {})
location = place.get('location')
if location:
lon, lat = location
if is_in_masovian(lat, lon):
city = place.get('city', 'Unknown')
if city in excluded_cities:
# Skip this station
continue
for module, module_data in station.get('measures', {}).items():
if 'temperature' in module_data.get('type', []):
res = module_data.get('res', {})
if res:
latest_timestamp = max(res.keys())
temperature = res[latest_timestamp][0]
stations.append(city)
temperatures.append(temperature)
coordinates.append((lat, lon))
break
return stations, temperatures, coordinates
def fetch_imgw_data():
url = "https://rafalraczynski.com.pl/imgw/dane-imgw/getJSON.php?type=table&province=14&sort=temp&order=asc"
response = requests.get(url)
if response.status_code == 200:
data = response.json()
else:
return [], [], []
target_stations = [
"ANDRZEJEWO", "BORKOWO", "CZARNOWO", "KAZANÓW",
"RUSZKOWICE", "RYBIENKO", "RZĄŚNIK WŁOŚCIAŃSKI", "WIELGOLAS"
]
stations = []
temperatures = []
station_ids = []
for item in data:
statName = item.get('statName')
temp = item.get('temp')
if statName and temp and statName.upper() in target_stations:
stations.append(statName)
temperatures.append(temp)
station_ids.append(item.get('statId', 'Unknown'))
return stations, temperatures, station_ids
def get_coordinates(station_name, geolocator):
try:
location = geolocator.geocode(f"{station_name}, Masovian Voivodeship, Poland")
if location:
if is_in_masovian(location.latitude, location.longitude):
return (location.latitude, location.longitude), False
else:
return None, False
return None, False
except Exception as e:
return None, False
def clean_temperature(temp):
cleaned_temp = re.sub(r'[↓↑\s]', '', temp)
try:
return float(cleaned_temp)
except ValueError:
return None
def remove_nearby_duplicates(trax_coords, trax_temperatures, netatmo_coords, netatmo_temperatures, distance_threshold=0.01):
trax_points = np.array([(lon, lat) for station, (lat, lon) in trax_coords.items()])
netatmo_points = np.array([(lon, lat) for lat, lon in netatmo_coords])
trax_tree = cKDTree(trax_points)
netatmo_filtered_coords = []
netatmo_filtered_temperatures = []
for netatmo_point, netatmo_temp in zip(netatmo_points, netatmo_temperatures):
distance, index = trax_tree.query(netatmo_point)
if distance > distance_threshold:
netatmo_filtered_coords.append(tuple(netatmo_point))
netatmo_filtered_temperatures.append(netatmo_temp)
return netatmo_filtered_coords, netatmo_filtered_temperatures
def hillshade(array, azimuth=315, angle_altitude=45):
# Convert angles to radians
azimuth_rad = np.radians(azimuth)
altitude_rad = np.radians(angle_altitude)
# Calculate gradients in x and y directions
x, y = np.gradient(array)
# Calculate the slope and aspect
slope = np.pi/2.0 - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
# Calculate shaded relief
shaded = np.sin(altitude_rad) * np.sin(slope) + np.cos(altitude_rad) * np.cos(slope) * np.cos(azimuth_rad - aspect)
# Normalize the shaded relief
shaded = (shaded + 1) / 2
shaded = (shaded * 255).astype(np.uint8)
return shaded
def read_and_process_modis_data(modis_day_folder, modis_night_folder):
# Read and stack MODIS day data
modis_day_files = sorted(glob.glob(os.path.join(modis_day_folder, '*.tif')))
modis_day_arrays = []
modis_day_transform = None
for modis_file in modis_day_files:
with rasterio.open(modis_file) as src:
data = src.read(1)
# Replace fill values with np.nan
fill_value = src.nodata
if fill_value is None:
# If nodata is not set, check for common fill values
data = np.where((data == -9999) | (data <= -1e20), np.nan, data)
else:
data = np.where(data == fill_value, np.nan, data)
modis_day_arrays.append(data)
if modis_day_transform is None:
modis_day_transform = src.transform
if modis_day_arrays:
modis_day_stack = np.stack(modis_day_arrays)
modis_day_data = np.nanmean(modis_day_stack, axis=0)
else:
modis_day_data = None
# Read and stack MODIS night data
modis_night_files = sorted(glob.glob(os.path.join(modis_night_folder, '*.tif')))
modis_night_arrays = []
modis_night_transform = None
for modis_file in modis_night_files:
with rasterio.open(modis_file) as src:
data = src.read(1)
# Replace fill values with np.nan
fill_value = src.nodata
if fill_value is None:
data = np.where((data == -9999) | (data <= -1e20), np.nan, data)
else:
data = np.where(data == fill_value, np.nan, data)
modis_night_arrays.append(data)
if modis_night_transform is None:
modis_night_transform = src.transform
if modis_night_arrays:
modis_night_stack = np.stack(modis_night_arrays)
modis_night_data = np.nanmean(modis_night_stack, axis=0)
else:
modis_night_data = None
return (modis_day_data, modis_day_transform), (modis_night_data, modis_night_transform)
def resample_modis_data(modis_data, modis_transform, dest_shape, dest_transform):
if modis_data is None:
return None
else:
# Resample to match the destination grid
dest_modis_data = np.zeros(dest_shape, dtype=modis_data.dtype)
reproject(
source=modis_data,
destination=dest_modis_data,
src_transform=modis_transform,
src_crs='EPSG:4326',
dst_transform=dest_transform,
dst_crs='EPSG:4326',
resampling=Resampling.nearest # Use nearest to maintain original resolution
)
# Replace zeros (from areas without data) with NaN
dest_modis_data = np.where(dest_modis_data == 0, np.nan, dest_modis_data)
return dest_modis_data
def plot_stations_with_boundary_and_interpolation(
trax_stations, trax_temperatures, trax_coords,
netatmo_stations, netatmo_temperatures, netatmo_coords,
imgw_stations, imgw_temperatures, imgw_coords,
modis_day_folder, modis_night_folder):
# Determine the current time in CET
cet = pytz.timezone('CET')
now_cet = datetime.now(cet)
hour = now_cet.hour + now_cet.minute / 60.0
# Define daytime and nighttime hours
day_start = 6
day_end = 18
is_daytime = day_start <= hour <= day_end
# List of Traxelektronik stations to exclude during daytime
exclude_trax_stations_daytime = [
'Ceranów',
'Baranów',
'Drobin',
'Kamion',
'Mińsk Mazowiecki',
'Mostówka',
'Paplin',
'Pilawa',
'Stromiec',
'Węzeł Łabiszyńska',
'Węzeł Marki',
'Węzeł Marynarska',
'Wilchta',
'Wiśniew'
]
# Clean the exclude list for matching
exclude_trax_stations_daytime_cleaned = [clean_station_name(name) for name in exclude_trax_stations_daytime]
# Map station names to cleaned station names
cleaned_station_dict = {station: clean_station_name(station) for station in trax_stations}
# Process Traxelektronik data
trax_filtered_stations = []
trax_filtered_temperatures = []
trax_geometries = []
# Filter and clean Traxelektronik data
for station, temperature in zip(trax_stations, trax_temperatures):
cleaned_station = cleaned_station_dict[station]
# Exclude certain stations during daytime
if is_daytime and cleaned_station in exclude_trax_stations_daytime_cleaned:
continue # Skip this station
cleaned_temperature = clean_temperature(temperature)
if cleaned_station in trax_coords and cleaned_temperature is not None:
lat, lon = trax_coords[cleaned_station]
trax_filtered_stations.append(cleaned_station)
trax_filtered_temperatures.append(cleaned_temperature)
trax_geometries.append(Point(lon, lat))
# Create GeoDataFrame for Traxelektronik data
trax_gdf = gpd.GeoDataFrame({
'Station': trax_filtered_stations,
'Temperature': trax_filtered_temperatures,
'Source': 'Traxelektronik'
}, geometry=trax_geometries, crs="EPSG:4326")
# Handle Netatmo data
netatmo_filtered_coords = []
netatmo_filtered_temperatures = []
netatmo_filtered_stations = []
for station, temperature, coord in zip(netatmo_stations, netatmo_temperatures, netatmo_coords):
if isinstance(temperature, (int, float)) and coord is not None:
netatmo_filtered_coords.append(coord)
netatmo_filtered_temperatures.append(temperature)
netatmo_filtered_stations.append(station)
if len(netatmo_filtered_coords) > 0:
netatmo_geometries = [Point(lon, lat) for lat, lon in netatmo_filtered_coords]
netatmo_gdf = gpd.GeoDataFrame({
'Station': netatmo_filtered_stations,
'Temperature': netatmo_filtered_temperatures,
'Source': 'Netatmo'
}, geometry=netatmo_geometries, crs="EPSG:4326")
else:
netatmo_gdf = gpd.GeoDataFrame(columns=['Station', 'Temperature', 'Source', 'geometry'], crs="EPSG:4326")
# Handle IMGW data
imgw_filtered_stations = []
imgw_filtered_temperatures = []
imgw_geometries = []
# Clean and process IMGW data
for station, temperature in zip(imgw_stations, imgw_temperatures):
cleaned_temperature = clean_temperature(temperature)
if station in imgw_coords and cleaned_temperature is not None:
lat, lon = imgw_coords[station]
imgw_filtered_stations.append(station)
imgw_filtered_temperatures.append(cleaned_temperature)
imgw_geometries.append(Point(lon, lat))
# Create GeoDataFrame for IMGW data
imgw_gdf = gpd.GeoDataFrame({
'Station': imgw_filtered_stations,
'Temperature': imgw_filtered_temperatures,
'Source': 'IMGW'
}, geometry=imgw_geometries, crs="EPSG:4326")
# Add X and Y coordinates to imgw_gdf
imgw_gdf['X'] = imgw_gdf.geometry.x
imgw_gdf['Y'] = imgw_gdf.geometry.y
# Read and clip DEM
dem_data, dem_transform = read_and_clip_dem(masovian_boundary)
# Read and clip Forest Cover
forest_data, forest_transform = read_and_clip_forest(masovian_boundary)
# Read and clip Water Bodies
water_data, water_transform, water_crs = read_and_clip_water_bodies(masovian_boundary)
# Ensure geometries are in EPSG:4326
trax_gdf_dem_crs = trax_gdf.to_crs("EPSG:4326")
netatmo_gdf_dem_crs = netatmo_gdf.to_crs("EPSG:4326")
imgw_gdf_dem_crs = imgw_gdf.to_crs("EPSG:4326")
# Get elevation at station points using dem_data and dem_transform
trax_elevations = get_elevation_at_points(
trax_gdf_dem_crs.geometry,
dem_data=dem_data,
dem_transform=dem_transform
)
netatmo_elevations = get_elevation_at_points(
netatmo_gdf_dem_crs.geometry,
dem_data=dem_data,
dem_transform=dem_transform
)
imgw_elevations = get_elevation_at_points(
imgw_gdf_dem_crs.geometry,
dem_data=dem_data,
dem_transform=dem_transform
)
# Add elevation to the GeoDataFrames
trax_gdf['Elevation'] = trax_elevations
netatmo_gdf['Elevation'] = netatmo_elevations
imgw_gdf['Elevation'] = imgw_elevations
# Get tree cover density at station points using forest_data and forest_transform
trax_forest = get_raster_value_at_points(
trax_gdf_dem_crs.geometry,
raster_data=forest_data,
raster_transform=forest_transform
)
netatmo_forest = get_raster_value_at_points(
netatmo_gdf_dem_crs.geometry,
raster_data=forest_data,
raster_transform=forest_transform
)
imgw_forest = get_raster_value_at_points(
imgw_gdf_dem_crs.geometry,
raster_data=forest_data,
raster_transform=forest_transform
)
# Add tree cover density to the GeoDataFrames
trax_gdf['Tree_Cover'] = trax_forest
netatmo_gdf['Tree_Cover'] = netatmo_forest
imgw_gdf['Tree_Cover'] = imgw_forest
# Get water body data at station points
trax_water = get_water_body_at_points(
trax_gdf_dem_crs.geometry,
water_data=water_data,
water_transform=water_transform
)
netatmo_water = get_water_body_at_points(
netatmo_gdf_dem_crs.geometry,
water_data=water_data,
water_transform=water_transform
)
imgw_water = get_water_body_at_points(
imgw_gdf_dem_crs.geometry,
water_data=water_data,
water_transform=water_transform
)
# Add water body data to the GeoDataFrames
trax_gdf['Water_Body'] = trax_water
netatmo_gdf['Water_Body'] = netatmo_water
imgw_gdf['Water_Body'] = imgw_water
# Add 'Priority' column to each GeoDataFrame
trax_gdf['Priority'] = 2
netatmo_gdf['Priority'] = 1
imgw_gdf['Priority'] = 3 # Higher priority
# Combine GeoDataFrames
stations_gdf = gpd.GeoDataFrame(pd.concat([trax_gdf, netatmo_gdf, imgw_gdf], ignore_index=True))
stations_gdf.crs = "EPSG:4326"
stations_gdf = stations_gdf.dropna(subset=['Temperature', 'Elevation', 'Tree_Cover'])
# Extract coordinates
stations_gdf['X'] = stations_gdf.geometry.x
stations_gdf['Y'] = stations_gdf.geometry.y
# Read ECOSTRESS data
ecostress_data, ecostress_transform = read_and_clip_ecostress(masovian_boundary)
# Read MODIS data
(modis_day_data, modis_day_transform), (modis_night_data, modis_night_transform) = read_and_process_modis_data(
modis_day_folder, modis_night_folder)
# Adjust grid_spacing as per your requirements
grid_spacing = 0.0025 # Approximately ~250 meters
x_min, y_min, x_max, y_max = masovian_boundary.total_bounds
# Compute the number of cols and rows with the new grid_spacing
ncols = int(np.ceil((x_max - x_min) / grid_spacing))
nrows = int(np.ceil((y_max - y_min) / grid_spacing))
# Create the affine transform for the grid
transform = from_origin(x_min, y_max, grid_spacing, grid_spacing)
# Create arrays of row and column indices
cols = np.arange(ncols)
rows = np.arange(nrows)
# Create a meshgrid of indices
cols_indices, rows_indices = np.meshgrid(cols, rows)
# Flatten the indices
cols_flat = cols_indices.ravel()
rows_flat = rows_indices.ravel()
# Get the x and y coordinates for each grid cell center
xs_flat, ys_flat = rasterio.transform.xy(transform, rows_flat, cols_flat)
# Resample MODIS data to match grid
modis_day_resampled = resample_modis_data(
modis_day_data, modis_day_transform,
dest_shape=(nrows, ncols), dest_transform=transform
)
modis_night_resampled = resample_modis_data(
modis_night_data, modis_night_transform,
dest_shape=(nrows, ncols), dest_transform=transform
)
# Decide which MODIS data to use based on time of day
if is_daytime:
modis_resampled = modis_day_resampled
else:
modis_resampled = modis_night_resampled
# After resampling, replace any extreme negative values with np.nan
modis_resampled = np.where(modis_resampled <= -1e20, np.nan, modis_resampled)
# Include MODIS LST data at station points
stations_gdf['MODIS_LST'] = get_raster_value_at_points(
stations_gdf.geometry,
raster_data=modis_resampled,
raster_transform=transform
)
# Handle missing MODIS_LST values
# Compute mean over finite values
modis_lst_mean = np.nanmean(stations_gdf['MODIS_LST'])
stations_gdf['MODIS_LST'].fillna(modis_lst_mean, inplace=True)
# Add MODIS_LST as a covariate
X = stations_gdf[['X', 'Y', 'Elevation', 'Tree_Cover', 'Water_Body', 'MODIS_LST']]
# Replace infinite values with NaN
X = X.replace([np.inf, -np.inf], np.nan)
# Remove rows with NaN values in X
X = X.dropna()
# Update y accordingly
y = stations_gdf.loc[X.index, 'Temperature']
# Proceed with scaling and model fitting
# Training Model Progress Bar
with tqdm(total=100, desc='Training Model', unit='%', ncols=80) as pbar:
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
pbar.update(30)
rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
rf_model.fit(X_scaled, y)
pbar.update(60)
# Get feature importances
feature_importances = rf_model.feature_importances_
feature_names = X.columns
# Create a DataFrame for better visualization
feature_importances_df = pd.DataFrame({'Feature': feature_names, 'Importance': feature_importances})
feature_importances_df = feature_importances_df.sort_values('Importance', ascending=False)
pbar.update(10)
# Create grid dataframe
grid_df = pd.DataFrame({
'X': xs_flat,
'Y': ys_flat,
'row_idx': rows_flat,
'col_idx': cols_flat
})
# Running Model Progress Bar
with tqdm(total=100, desc='Running Model', unit='%', ncols=80) as pbar:
# Retrieve elevations from DEM
grid_elevations = get_elevation_at_points(
[Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
dem_data=dem_data,
dem_transform=dem_transform
)
grid_df['Elevation'] = grid_elevations
pbar.update(20)
# Retrieve tree cover density from forest cover raster
grid_forest = get_raster_value_at_points(
[Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
raster_data=forest_data,
raster_transform=forest_transform
)
grid_df['Tree_Cover'] = grid_forest
pbar.update(20)
# Retrieve water body data at grid points
grid_water = get_water_body_at_points(
[Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
water_data=water_data,
water_transform=water_transform
)
grid_df['Water_Body'] = grid_water
pbar.update(20)
# Remove grid points with NaN values (i.e., outside data coverage)
grid_df = grid_df.dropna(subset=['Elevation', 'Tree_Cover'])
# Include MODIS_LST in grid_df
grid_df['MODIS_LST'] = grid_df.apply(
lambda row: get_raster_value_at_point(row['X'], row['Y'], modis_resampled, transform),
axis=1
)
# Handle missing MODIS_LST values
modis_lst_mean_grid = np.nanmean(grid_df['MODIS_LST'])
grid_df['MODIS_LST'].fillna(modis_lst_mean_grid, inplace=True)
pbar.update(20)
# Prepare data for regression prediction
X_grid = grid_df[['X', 'Y', 'Elevation', 'Tree_Cover', 'Water_Body', 'MODIS_LST']]
# Handle infinite values in X_grid
X_grid = X_grid.replace([np.inf, -np.inf], np.nan)
# Handle NaNs in X_grid
X_grid = X_grid.fillna(method='ffill').fillna(method='bfill')
# Ensure that all values are finite
# Scale X_grid
X_grid_scaled = scaler.transform(X_grid)
pbar.update(10)
grid_df['Predicted_Temperature'] = rf_model.predict(X_grid_scaled)
pbar.update(10)
# Calculate residuals at station locations
stations_gdf['Predicted_Temperature'] = rf_model.predict(X_scaled)
stations_gdf['Residual'] = stations_gdf['Temperature'] - stations_gdf['Predicted_Temperature']
# Adjust anomalously warm values from weather stations
# Calculate percentiles of residuals
lower_percentile = np.percentile(stations_gdf['Residual'], 2.5)
upper_percentile = np.percentile(stations_gdf['Residual'], 97.5)
# Cap the residuals
stations_gdf['Adjusted_Residual'] = stations_gdf['Residual'].clip(lower=lower_percentile, upper=upper_percentile)
# Remove duplicates by prioritizing IMGW data
data = stations_gdf.copy()
# Sort by priority (higher priority first)
data = data.sort_values('Priority', ascending=False)
# Group by coordinates and keep the first record in each group (highest priority)
data = data.groupby(['X', 'Y'], as_index=False).first()
# Set up variogram model with adjusted parameters
variogram_model = 'spherical'
variogram_parameters = {'nugget': 0.2, 'sill': 1.0, 'range': 0.05}
# Perform Universal Kriging using Adjusted_Residual
uk = UniversalKriging(
x=data['X'].values,
y=data['Y'].values,
z=data['Adjusted_Residual'].values,
variogram_model=variogram_model,
variogram_parameters=variogram_parameters,
nlags=20,
verbose=False
)
z, ss = uk.execute(
'points',
grid_df['X'].values,
grid_df['Y'].values
)
grid_df['Residual'] = z
grid_df['Kriging_Variance'] = ss
# Handle NaNs in residuals by filling them with zero
grid_df['Residual'].fillna(0, inplace=True)
# Final temperature estimation
grid_df['Temperature'] = grid_df['Predicted_Temperature'] + grid_df['Residual']
# Save a copy of the temperature before ECOSTRESS adjustments
grid_df['Temperature_No_ECOSTRESS'] = grid_df['Temperature'].copy()
# Topographic Cooling Effect
max_cooling_effect = 3.0 # Maximum cooling effect in degrees Celsius
# Normalize elevation to range from 0 to 1
elevation_min = grid_df['Elevation'].min()
elevation_max = grid_df['Elevation'].max()
grid_df['Normalized_Elevation'] = (grid_df['Elevation'] - elevation_min) / (elevation_max - elevation_min)
# Compute the topographic cooling effect
grid_df['Topographic_Cooling'] = grid_df['Normalized_Elevation'] * max_cooling_effect * (1 - is_daytime)
# Ensure topographic cooling is not applied over water bodies
grid_df.loc[grid_df['Water_Body'] == 1, 'Topographic_Cooling'] = 0
# Apply the topographic cooling effect to the temperature
grid_df['Temperature'] -= grid_df['Topographic_Cooling']
# ECOSTRESS Adjustments
# Retrieve ECOSTRESS LST values at grid points (same as before)
grid_ecostress = get_raster_value_at_points(
[Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
raster_data=ecostress_data,
raster_transform=ecostress_transform
)
grid_df['ECOSTRESS_LST'] = grid_ecostress
# Handle NaNs in ECOSTRESS_LST by filling with MODIS_LST values or the mean value
grid_df['ECOSTRESS_LST'] = grid_df.apply(
lambda row: row['ECOSTRESS_LST'] if not np.isnan(row['ECOSTRESS_LST']) else row['MODIS_LST'], axis=1
)
grid_df['ECOSTRESS_LST'].fillna(grid_df['ECOSTRESS_LST'].mean(), inplace=True)
# Compute ECOSTRESS anomalies (deviation from mean)
ecostress_mean = grid_df['ECOSTRESS_LST'].mean()
grid_df['ECOSTRESS_Anomaly'] = grid_df['ECOSTRESS_LST'] - ecostress_mean
# --- Begin Anomaly Adjustment for Nighttime Forested Areas ---
if not is_daytime:
forest_threshold = 50 # Define the threshold for forest cover (e.g., 50%)
warm_anomaly_forest_mask = (grid_df['ECOSTRESS_Anomaly'] > 0) & (grid_df['Tree_Cover'] >= forest_threshold)
# Invert the anomalies in forested areas with warm anomalies at night
grid_df.loc[warm_anomaly_forest_mask, 'ECOSTRESS_Anomaly'] *= -1
# --- End Anomaly Adjustment ---
# Normalize anomalies to range [-1, 1]
ecostress_anomaly_min = grid_df['ECOSTRESS_Anomaly'].min()
ecostress_anomaly_max = grid_df['ECOSTRESS_Anomaly'].max()
grid_df['ECOSTRESS_Anomaly_Normalized'] = (
(grid_df['ECOSTRESS_Anomaly'] - ecostress_anomaly_min) / (ecostress_anomaly_max - ecostress_anomaly_min) * 2 - 1
)
# Increase the impact of ECOSTRESS LST data
amplification_factor = 3.0 # Adjust as needed
# Compute night_weight: maximum at midnight, zero at noon
night_weight = (math.cos(((hour % 24) / 24) * 2 * math.pi) + 1) / 2 # Ranges from 1 at midnight to 0 at noon
# Apply anomalies with the night_weight
adjustment = grid_df['ECOSTRESS_Anomaly_Normalized'] * amplification_factor * night_weight
# Prevent ECOSTRESS adjustments on water bodies
adjustment[grid_df['Water_Body'] == 1] = 0
# Limit maximum adjustment to avoid over-adjusting
max_adjustment = 5.0 # degrees Celsius
adjustment = adjustment.clip(-max_adjustment, max_adjustment)
# Apply adjustments to the interpolated temperatures
grid_df['Temperature'] += adjustment
# Apply cooling effect on water bodies
cooling_effect = -2.0 # degrees Celsius
grid_df.loc[grid_df['Water_Body'] == 1, 'Temperature'] += cooling_effect
# Apply Gaussian smoothing to the interpolated temperature grid
# Reintroduce smoothing to reduce square artifacts
from scipy.ndimage import gaussian_filter
# Create 2D arrays of temperatures
grid_temperature_with_ecostress = np.full((nrows, ncols), np.nan)
grid_temperature_with_ecostress[grid_df['row_idx'], grid_df['col_idx']] = grid_df['Temperature'].values
# Apply Gaussian smoothing
sigma = 1.0 # Adjust sigma to control the amount of smoothing
temperature_with_ecostress_smoothed = gaussian_filter(grid_temperature_with_ecostress, sigma=sigma)
# Mask the grid_temperature with the Masovian boundary
masovian_shape = [masovian_boundary.geometry.unary_union]
mask = rasterize(
masovian_shape,
out_shape=(nrows, ncols),
transform=transform,
fill=0,
all_touched=True,
default_value=1,
dtype='uint8'
)
# Apply the mask
temperature_with_ecostress_smoothed = np.ma.array(temperature_with_ecostress_smoothed, mask=mask == 0)
# Use the smoothed temperature grid
temperature_with_ecostress = temperature_with_ecostress_smoothed
# Correct anomalously warm interpolated values based on nearest IMGW stations
# Build a KDTree for IMGW station locations
imgw_coords_array = imgw_gdf[['X', 'Y']].values
imgw_temperatures_array = imgw_gdf['Temperature'].values
imgw_tree = cKDTree(imgw_coords_array)
# For each grid point, find the nearest IMGW station
grid_points = grid_df[['X', 'Y']].values
distances, indices = imgw_tree.query(grid_points, k=1)
# Get the temperatures of the nearest IMGW stations
nearest_imgw_temperatures = imgw_temperatures_array[indices]
# Identify grid points where the interpolated temperature is significantly higher than the nearest IMGW station
temp_difference = temperature_with_ecostress.flatten() - nearest_imgw_temperatures
anomaly_threshold = 2.0 # degrees Celsius
anomalous_warm_mask = temp_difference > anomaly_threshold
# Apply correction to anomalously warm grid points
# Reduce the temperature difference to the threshold
corrected_temperatures = temperature_with_ecostress.flatten()
corrected_temperatures[anomalous_warm_mask] = nearest_imgw_temperatures[anomalous_warm_mask] + anomaly_threshold
# Reshape corrected_temperatures back to the grid shape
corrected_temperature_grid = corrected_temperatures.reshape((nrows, ncols))
corrected_temperature_grid = np.ma.array(corrected_temperature_grid, mask=mask == 0)
# Update the temperature_with_ecostress with corrected temperatures
temperature_with_ecostress = corrected_temperature_grid
# Also create grid_variance array (we'll skip smoothing variance for simplicity)
grid_variance = np.full((nrows, ncols), np.nan)
grid_variance[grid_df['row_idx'], grid_df['col_idx']] = grid_df['Kriging_Variance']
grid_variance = np.ma.array(grid_variance, mask=mask == 0)
# Compute Standard Deviation from Variance
grid_std_dev = np.sqrt(grid_variance)
grid_std_dev = np.ma.array(grid_std_dev, mask=mask == 0)
# Read the custom color scale from the CSV file
color_scale_df = pd.read_csv('input/color_scale.csv')
color_scale_df = color_scale_df.sort_values('value')
values = color_scale_df['value'].values
hex_colors = color_scale_df['color'].values
# Ensure colors start with '#'
hex_colors = [color if color.startswith('#') else f"#{color}" for color in hex_colors]
# Create a custom colormap
# Normalize the values to 0..1
vmin = values.min()
vmax = values.max()
normalized_values = (values - vmin) / (vmax - vmin)
# Create a list of tuples (position, color)
color_tuples = list(zip(normalized_values, hex_colors))
# Create a LinearSegmentedColormap
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', color_tuples)
# Set the normalization
norm = mcolors.Normalize(vmin=vmin, vmax=vmax)
with rasterio.open('input/masovian_dem.tif') as dem_src:
# Reproject Masovian boundary to DEM CRS if necessary
if masovian_boundary.crs != dem_src.crs:
masovian_boundary_dem_crs = masovian_boundary.to_crs(dem_src.crs)
else:
masovian_boundary_dem_crs = masovian_boundary
# Clip the DEM to the Masovian boundary
dem_plot_data, dem_plot_transform = rasterio.mask.mask(
dem_src, masovian_boundary_dem_crs.geometry, crop=True, filled=True, nodata=np.nan
)
# Read the first band (assuming single-band DEM)
dem_plot_data = dem_plot_data[0] # Since mask returns data as (bands, rows, cols)
# Handle NaNs in the DEM data
dem_plot_data = np.nan_to_num(dem_plot_data, nan=np.nanmin(dem_plot_data))
# Calculate the extent of the clipped DEM for plotting
dem_left, dem_top = dem_plot_transform * (0, 0)
dem_right, dem_bottom = dem_plot_transform * (dem_plot_data.shape[1], dem_plot_data.shape[0])
dem_extent = [dem_left, dem_right, dem_bottom, dem_top]
# Compute extent for plotting (same as before)
def get_extent(transform, array):
x_min = transform.c
y_max = transform.f
x_max = x_min + transform.a * array.shape[1]
y_min = y_max + transform.e * array.shape[0]
return [x_min, x_max, y_min, y_max]
hs = hillshade(dem_plot_data, azimuth=315, angle_altitude=45)
extent = get_extent(transform, temperature_with_ecostress)
dem_extent = get_extent(dem_plot_transform, dem_plot_data)
# Read water data and get CRS
water_data, water_transform, water_crs = read_and_clip_water_bodies(masovian_boundary)
# Create an empty array to hold the resampled water data
water_resampled = np.empty((nrows, ncols), dtype=water_data.dtype)
# Perform the reprojection and resampling
reproject(
source=water_data,
destination=water_resampled,
src_transform=water_transform,
src_crs=water_crs,
dst_transform=transform,
dst_crs='EPSG:4326', # Assuming temperature grid is in EPSG:4326
resampling=Resampling.nearest
)
# Now water_resampled is aligned with the temperature grid
# Create the water mask
water_mask = np.where(water_resampled > 0, 1, 0)
# Define the bounding box coordinates (EPSG:4326)
xmin, ymin, xmax, ymax = 20.188477, 51.017203, 22.061646, 51.821759
# Compute the indices for the bounding box
from rasterio.transform import rowcol
# Convert coordinates to row and column indices in the temperature grid
row_ul, col_ul = rowcol(transform, xmin, ymax)
row_lr, col_lr = rowcol(transform, xmax, ymin)
# Ensure indices are within valid ranges
nrows, ncols = temperature_with_ecostress.shape
row_ul = max(0, min(nrows - 1, row_ul))
row_lr = max(0, min(nrows - 1, row_lr))
col_ul = max(0, min(ncols - 1, col_ul))
col_lr = max(0, min(ncols - 1, col_lr))
# Ensure row_ul <= row_lr and col_ul <= col_lr
if row_ul > row_lr:
row_ul, row_lr = row_lr, row_ul
if col_ul > col_lr:
col_ul, col_lr = col_lr, col_ul
# Extract the subset arrays
temperature_subset = temperature_with_ecostress[row_ul:row_lr+1, col_ul:col_lr+1]
hs_subset = hs[row_ul:row_lr+1, col_ul:col_lr+1]
water_mask_subset = water_mask[row_ul:row_lr+1, col_ul:col_lr+1]
mask_subset = mask[row_ul:row_lr+1, col_ul:col_lr+1]
# Apply the Masovian mask to the subsets
temperature_subset = np.ma.array(temperature_subset, mask=mask_subset == 0)
hs_subset = np.ma.array(hs_subset, mask=mask_subset == 0)
water_mask_subset = np.ma.array(water_mask_subset, mask=mask_subset == 0)
grid_spacing_x = transform.a # Positive grid spacing in x direction
grid_spacing_y = transform.e # Negative grid spacing in y direction
new_x_min = transform.c + col_ul * grid_spacing_x
new_y_max = transform.f + row_ul * grid_spacing_y # grid_spacing_y is negative
new_transform = from_origin(new_x_min, new_y_max, grid_spacing_x, -grid_spacing_y)
# Compute the new extent
extent_subset = get_extent(new_transform, temperature_subset)
# Filter stations within the bounding box
from shapely.geometry import box
bbox_polygon = box(xmin, ymin, xmax, ymax)
stations_in_bbox = stations_gdf[stations_gdf.geometry.within(bbox_polygon)]
# Clip the Masovian boundary to the bounding box
masovian_boundary_clipped = masovian_boundary.clip(bbox_polygon)
# Plotting Temperature Map with Contours
# Calculate the aspect ratio of the bounding box
width = xmax - xmin
height = ymax - ymin
aspect_ratio = width / height
# Set the figure size based on the aspect ratio
desired_height_inches = 10 # You can adjust this value as desired
desired_width_inches = desired_height_inches * aspect_ratio
# Plotting Temperature Map with ECOSTRESS Adjustments
fig1, ax1 = plt.subplots(figsize=(15, 15))
# Plot the interpolated temperature grid with ECOSTRESS adjustments
img1 = ax1.imshow(temperature_with_ecostress, extent=extent, origin='upper',
cmap=cmap, norm=norm, alpha=1)
# Plot hillshade on top
ax1.imshow(hs, extent=dem_extent, cmap='gray', origin='upper', alpha=0.3)
ax1.imshow(water_mask, extent=extent, origin='upper', cmap='Blues', alpha=0.5) # Increased alpha to 0.5
# Plot Masovian Voivodeship boundary with a light, dark stroke
masovian_boundary.boundary.plot(ax=ax1, edgecolor='#333333', linewidth=1.5)
# Plot stations
stations_gdf.plot(ax=ax1, color='black', marker='s', markersize=5)
# Remove axes for a minimalistic look
ax1.axis('off')
# Draw a light, dark stroke around the plot extent
from matplotlib.patches import Rectangle
rect = Rectangle(
(extent[0], extent[2]), # Lower-left corner (x0, y0)
extent[1] - extent[0], # Width
extent[3] - extent[2], # Height
linewidth=1.5,
edgecolor='#333333',
facecolor='none',
transform=ax1.transData
)
ax1.add_patch(rect)
# Introduce a subtle background gradient
from matplotlib.colors import LinearSegmentedColormap
background_cmap = LinearSegmentedColormap.from_list(
'background_gradient', ['#ffffff', '#f0f0f0'], N=256)
ax1.imshow(
np.linspace(0, 1, 256).reshape(1, -1),
cmap=background_cmap,
extent=(extent[0], extent[1], extent[2], extent[3]),
aspect='auto',
alpha=0.2,
zorder=-1 # Place it below other elements
)
# Adjust font properties for a modern look
plt.rcParams.update({
'font.size': 14, # Base font size
'font.family': 'sans-serif',
'font.sans-serif': ['Arial', 'Helvetica', 'DejaVu Sans'],
'text.color': '#333333'
})
# Enlarge and reposition the product name (Title)
fig1.text(
0.5, 0.94, # Adjusted y-position for better placement
"High-Resolution Mesoscale Temperature Analysis",
fontsize=28, # Increased font size
fontweight='bold',
ha='center',
va='top',
fontname='DejaVu Sans',
color='#333333'
)
# Move the model name/resolution text closer to the product name text
fig1.text(
0.5, 0.91, # Adjusted y-position to move it closer
"Model: TempMeso v1.0 (Preview) | Resolution: 250 m",
fontsize=18, # Increased font size
ha='center',
va='top',
fontname='DejaVu Sans',
color='#666666'
)
# Adjust colorbar position and style
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="3%", pad=0.1)
cbar1 = fig1.colorbar(img1, cax=cax)
ticks = np.arange(np.floor(vmin / 5) * 5, np.ceil(vmax / 5) * 5 + 1, 5)
cbar1.set_ticks(ticks)
cbar1.set_ticklabels(ticks.astype(int))
cbar1.ax.tick_params(labelsize=14)
# Enhance colorbar aesthetics
cbar1.outline.set_visible(False)
cbar1.ax.yaxis.set_tick_params(width=0.5)
cbar1.ax.tick_params(labelsize=14, colors='#333333')
# Move the logo to the bottom right corner, directly below the colorbar
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
# Read the logo image
logo_path = 'input/logo.png'
logo_img = mpimg.imread(logo_path)
# Adjust the size of the logo
logo_zoom = 0.25
imagebox = OffsetImage(logo_img, zoom=logo_zoom)
# Get the position of the colorbar axes
cbar_bbox = cax.get_position()
cbar_x0 = cbar_bbox.x0 # Left
cbar_x1 = cbar_bbox.x1 # Right
# Calculate the x position for the logo (centered below the colorbar)
logo_x = (cbar_x0 + cbar_x1) / 2
# Set the y position as low as possible, but ensure the logo is fully visible
logo_y = 0.02 # Adjust this value if necessary
ab = AnnotationBbox(
imagebox,
xy=(logo_x, logo_y),
xycoords='figure fraction',
frameon=False,
box_alignment=(0.5, 0), # Centered horizontally, aligned at bottom
)
fig1.add_artist(ab)
# Adjust the bottom margin to ensure the logo is not cut off
plt.subplots_adjust(left=0.05, right=0.95, top=0.88, bottom=0.05)
# On the map, add texts at the coldest and warmest temperature locations
# Calculate min and max temperatures over the entire Masovian region
min_temp = np.nanmin(temperature_with_ecostress)
max_temp = np.nanmax(temperature_with_ecostress)
# Get the indices of the min and max temperatures
min_indices = np.unravel_index(np.nanargmin(temperature_with_ecostress), temperature_with_ecostress.shape)
max_indices = np.unravel_index(np.nanargmax(temperature_with_ecostress), temperature_with_ecostress.shape)
# Convert grid indices to coordinates
min_lon_full, min_lat_full = rasterio.transform.xy(transform, min_indices[0], min_indices[1])
max_lon_full, max_lat_full = rasterio.transform.xy(transform, max_indices[0], max_indices[1])
# Plot the temperature values on the map (without markers)
# Use white text with a thick black stroke
text_kwargs = dict(
fontsize=18, # Increased font size
fontname='Arial',
ha='center',
va='center',
color='white',
transform=ax1.transData,
path_effects=[patheffects.withStroke(linewidth=3, foreground='black')]
)
# Format temperature values with one decimal point
ax1.text(min_lon_full, min_lat_full, f"{min_temp:.1f}", **text_kwargs)
ax1.text(max_lon_full, max_lat_full, f"{max_temp:.1f}", **text_kwargs)
# Save and show the plot
plt.savefig('output/temperature_analysis_modern.png', dpi=300)
plt.show()
print("Updated temperature analysis plot created successfully.")
# Calculate the aspect ratio of the bounding box
width = xmax - xmin
height = ymax - ymin
aspect_ratio = width / height
# Set the maximum desired width
max_width_inches = 15 # Adjust this value as desired
# Calculate the figure height based on the aspect ratio and maximum width
desired_width_inches = min(max_width_inches, aspect_ratio * desired_height_inches)
desired_height_inches = desired_width_inches / aspect_ratio
# Ensure the height doesn't exceed a certain value
desired_height_inches = min(desired_height_inches, 10) # Adjust as needed
fig2, ax2 = plt.subplots(figsize=(desired_width_inches, desired_height_inches))
# Proceed with previous plotting code, using adjusted fig2 and ax2
# Use the subset data and extent with applied mask
img2 = ax2.imshow(temperature_subset, extent=extent_subset, origin='upper',
cmap=cmap, norm=norm, alpha=1)
# Plot hillshade on top using subset data with mask
ax2.imshow(hs_subset, extent=extent_subset, cmap='gray', origin='upper', alpha=0.3)
# Overlay water bodies using subset data with mask
ax2.imshow(water_mask_subset, extent=extent_subset, origin='upper', cmap='Blues', alpha=0.5)
# Plot clipped Masovian boundary with a light, dark stroke
masovian_boundary_clipped.boundary.plot(ax=ax2, edgecolor='#333333', linewidth=1.5)
# Plot stations within the bounding box
stations_in_bbox.plot(ax=ax2, color='black', marker='s', markersize=5)
# Remove axes for a minimalistic look
ax2.axis('off')
# Introduce a subtle background gradient using the subset extent
ax2.imshow(
np.linspace(0, 1, 256).reshape(1, -1),
cmap=background_cmap,
extent=extent_subset,
aspect='auto',
alpha=0.2,
zorder=-1 # Place it below other elements
)
# Adjust font properties for a modern look
plt.rcParams.update({
'font.size': 14,
'font.family': 'sans-serif',
'font.sans-serif': ['Arial', 'Helvetica', 'DejaVu Sans'],
'text.color': '#333333'
})
# Enlarge and reposition the product name (Title)
fig2.text(
0.5, 0.995, # Moved up towards the top edge
"High-Resolution Mesoscale Temperature Analysis",
fontsize=12,
fontweight='bold',
ha='center',
va='top',
fontname='DejaVu Sans',
color='#333333'
)
# Move the model name/resolution text closer to the product name text
fig2.text(
0.5, 0.970, # Moved up slightly below the title
"Model: HRMTA v1.0 (Preview) | Resolution: 250 m",
fontsize=9,
ha='center',
va='top',
fontname='DejaVu Sans',
color='#666666'
)
# Add data source attribution at the bottom left corner
fig2.text(
0.01, 0.02,
"Data source: IMGW-PIB (processed)",
fontsize=12,
ha='left',
va='bottom',
fontname='DejaVu Sans',
color='#333333'
)
# Adjust colorbar position and style
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax2)
cax2 = divider.append_axes("right", size="3%", pad=0.1)
cbar2 = fig2.colorbar(img2, cax=cax2)
cbar2.set_ticks(ticks)
cbar2.set_ticklabels(ticks.astype(int))
# Import FontProperties and load the font
from matplotlib.font_manager import FontProperties
# Replace 'input/RobotoCondensed-Regular.ttf' with the actual path to your font file
roboto_font_path = 'input/RobotoCondensed-Regular.ttf'
roboto_condensed = FontProperties(fname=roboto_font_path)
# Apply font properties to the colorbar tick labels
for label in cbar2.ax.yaxis.get_ticklabels():
label.set_fontproperties(roboto_condensed)
cbar2.ax.tick_params(labelsize=14)
# Enhance colorbar aesthetics
cbar2.outline.set_visible(False)
cbar2.ax.yaxis.set_tick_params(width=0.5)
cbar2.ax.tick_params(labelsize=14, colors='#333333')
# Move the logo to the bottom right corner, directly below the colorbar
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
# Read the logo image again for the second plot
logo_img2 = mpimg.imread(logo_path)
# Adjust the size of the logo
imagebox2 = OffsetImage(logo_img2, zoom=logo_zoom)
# Get the position of the colorbar axes
cbar_bbox2 = cax2.get_position()
logo_x2 = (cbar_bbox2.x0 + cbar_bbox2.x1) / 2
logo_y2 = 0.02 # Adjust this value if necessary
ab2 = AnnotationBbox(
imagebox2,
xy=(logo_x2, logo_y2),
xycoords='figure fraction',
frameon=False,
box_alignment=(0.5, 0),
)
fig2.add_artist(ab2)
# Adjust the bottom margin to ensure the logo and data attribution text are not cut off
plt.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.06)
# Calculate min and max temperatures over the temperature_subset
min_temp = np.nanmin(temperature_subset)
max_temp = np.nanmax(temperature_subset)
# Get the indices of the min and max temperatures within the subset array
min_indices = np.unravel_index(np.nanargmin(temperature_subset), temperature_subset.shape)
max_indices = np.unravel_index(np.nanargmax(temperature_subset), temperature_subset.shape)
# Convert grid indices to coordinates using the new_transform
from rasterio.transform import xy
min_lon, min_lat = xy(new_transform, min_indices[0], min_indices[1])
max_lon, max_lat = xy(new_transform, max_indices[0], max_indices[1])
# Add contours to the plot using the subset data
min_temp_c = np.floor(min_temp)
max_temp_c = np.ceil(max_temp)
contour_levels = np.arange(min_temp_c, max_temp_c + 1, 1)
contours = ax2.contour(
temperature_subset, levels=contour_levels, colors='black', linewidths=0.5,
extent=extent_subset, origin='upper'
)
contour_labels = ax2.clabel(
contours, inline=False, fontsize=8, fmt='%d', colors='white',
inline_spacing=4, manual=False
)
for txt in contour_labels:
txt.set_path_effects([
patheffects.Stroke(linewidth=3, foreground='black'),
patheffects.Normal()
])
# Now add the min/max temperature labels after the contours and contour labels
text_kwargs = dict(
fontsize=18,
fontname='Arial',
ha='center',
va='center',
color='white',
transform=ax2.transData,
path_effects=[patheffects.withStroke(linewidth=3, foreground='black')]
)
# Format temperature values with one decimal point
ax2.text(min_lon, min_lat, f"{min_temp:.1f}", **text_kwargs, zorder=5)
ax2.text(max_lon, max_lat, f"{max_temp:.1f}", **text_kwargs, zorder=5)
# Reverse geocode min and max temperature locations to find nearest village
from geopy.geocoders import Nominatim
geolocator = Nominatim(user_agent="geoapiExercises")
try:
min_location = geolocator.reverse((min_lat_full, min_lon_full), language='en')
min_address = min_location.raw.get('address', {})
min_village = min_address.get('village',
min_address.get('town',
min_address.get('city', 'Unknown Location')))
except Exception as e:
print(f"Error reverse geocoding min location: {e}")
min_village = 'Unknown Location'
try:
max_location = geolocator.reverse((max_lat_full, max_lon_full), language='en')
max_address = max_location.raw.get('address', {})
max_village = max_address.get('village',
max_address.get('town',
max_address.get('city', 'Unknown Location')))
except Exception as e:
print(f"Error reverse geocoding max location: {e}")
max_village = 'Unknown Location'
# Create label texts
min_label_text = f"MIN\n {min_temp:.1f}°C: {min_village}"
max_label_text = f"MAX\n {max_temp:.1f}°C: {max_village}"
# Define font properties
from matplotlib.font_manager import FontProperties
roboto_font_path = 'input/RobotoCondensed-Regular.ttf'
roboto_condensed = FontProperties(fname=roboto_font_path)
# Add text labels to the figure
fig1.text(0.01, 0.95, min_label_text, ha='left', va='top', fontsize=14, fontproperties=roboto_condensed)
fig1.text(0.99, 0.95, max_label_text, ha='right', va='top', fontsize=14, fontproperties=roboto_condensed)
# Save and show the plot
plt.savefig('output/temperature_analysis_with_contours.png', dpi=300)
plt.show()
print("Temperature analysis plot with contours created successfully.")
# At the end, print feature importances to confirm MODIS influence
print("\nFeature importances from the Random Forest model:")
print(feature_importances_df)
# Check if MODIS_LST has significant importance
modis_importance = feature_importances_df.loc[feature_importances_df['Feature'] == 'MODIS_LST', 'Importance'].values[0]
if modis_importance > 0.05: # Threshold can be adjusted
print(f"\nMODIS_LST feature has significant influence on the model (importance: {modis_importance:.4f}).")
else:
print(f"\nMODIS_LST feature has low influence on the model (importance: {modis_importance:.4f}).")
return temperature_with_ecostress, transform
def get_elevation_at_points(geometry, dem_data=None, dem_transform=None):
from rasterio.transform import rowcol
if dem_data is not None and dem_transform is not None:
elevations = []
for point in geometry:
x, y = point.x, point.y
try:
row, col = rowcol(dem_transform, x, y)
row = int(row)
col = int(col)
elevation = dem_data[row, col]
elevations.append(elevation)
except (IndexError, ValueError):
elevations.append(np.nan)
return elevations
else:
raise ValueError("Provide dem_data and dem_transform.")
def get_interpolated_temperature_at_point(lat, lon, temperature_grid, transform):
"""
Given a latitude and longitude, return the interpolated temperature value at that point.
Parameters:
- lat, lon: Latitude and longitude of the point.
- temperature_grid: The interpolated temperature grid (2D numpy array).
- transform: Affine transform for the grid.
Returns:
- Interpolated temperature value at the point, or None if point is outside the grid.
"""
from rasterio.transform import rowcol
try:
row, col = rowcol(transform, lon, lat) # Note that rasterio uses (x, y), so (lon, lat)
row = int(row)
col = int(col)
value = temperature_grid[row, col]
if np.ma.is_masked(value) or np.isnan(value):
# Value is masked or NaN, meaning outside the interpolated area
return None
else:
return value.item()
except (IndexError, ValueError):
return None # Point is outside the grid
def get_raster_value_at_points(geometry, raster_data=None, raster_transform=None):
from rasterio.transform import rowcol
if raster_data is not None and raster_transform is not None:
values = []
for point in geometry:
x, y = point.x, point.y
try:
row, col = rowcol(raster_transform, x, y)
row = int(row)
col = int(col)
if (0 <= row < raster_data.shape[0]) and (0 <= col < raster_data.shape[1]):
value = raster_data[row, col]
# Check for NaN in raster data
if np.isnan(value):
values.append(np.nan)
else:
values.append(value)
else:
values.append(np.nan)
except (IndexError, ValueError):
values.append(np.nan)
return values
else:
raise ValueError("Provide raster_data and raster_transform.")
def get_raster_value_at_point(x, y, raster_data=None, raster_transform=None):
from rasterio.transform import rowcol
if raster_data is not None and raster_transform is not None:
try:
row, col = rowcol(raster_transform, x, y)
row = int(row)
col = int(col)
if (0 <= row < raster_data.shape[0]) and (0 <= col < raster_data.shape[1]):
value = raster_data[row, col]
if np.isnan(value):
return np.nan
else:
return value
else:
return np.nan
except (IndexError, ValueError):
return np.nan
else:
raise ValueError("Provide raster_data and raster_transform.")
def get_water_body_at_points(geometry, water_data=None, water_transform=None):
from rasterio.transform import rowcol
if water_data is not None and water_transform is not None:
water_bodies = []
for point in geometry:
x, y = point.x, point.y
try:
row, col = rowcol(water_transform, x, y)
row = int(row)
col = int(col)
if (0 <= row < water_data.shape[0]) and (0 <= col < water_data.shape[1]):
value = water_data[row, col]
# Adjust the condition based on your raster's water encoding
# For binary raster where 1 = water, 0 = land
water_bodies.append(1 if value else 0)
else:
water_bodies.append(0)
except (IndexError, ValueError):
water_bodies.append(0)
return water_bodies
else:
raise ValueError("Provide water_data and water_transform.")
def read_and_clip_dem(masovian_boundary):
dem_path = 'input/masovian_dem.tif'
with rasterio.open(dem_path) as src:
dem_data = src.read(1, masked=True)
dem_transform = src.transform
return dem_data, dem_transform
def read_and_clip_forest(masovian_boundary):
forest_path = 'input/forests.tif'
with rasterio.open(forest_path) as src:
forest_data = src.read(1, masked=True)
forest_data = np.where(forest_data == src.nodata, 0, forest_data)
forest_transform = src.transform
return forest_data, forest_transform
def read_and_clip_water_bodies(masovian_boundary):
water_path = 'input/water_bodies.tif'
with rasterio.open(water_path) as src:
water_crs = src.crs # Get the CRS of the water raster
# Ensure the Masovian boundary is in the same CRS as the water raster
if masovian_boundary.crs != src.crs:
masovian_boundary = masovian_boundary.to_crs(src.crs)
# Perform the clipping
clipped_water, water_transform = rasterio.mask.mask(
src, masovian_boundary.geometry, crop=True, filled=True, nodata=0
)
# Read the first band (assuming single-band raster)
water_data = clipped_water[0]
# Handle NoData values by setting them to 0 (assuming 0 represents no water)
water_data = np.where(water_data == src.nodata, 0, water_data)
return water_data, water_transform, water_crs
def read_and_clip_ecostress(masovian_boundary):
ecostress_path = 'input/raw_ecostress.tif'
with rasterio.open(ecostress_path) as src:
# Read data
ecostress_data = src.read(1, masked=True)
ecostress_data = np.where(ecostress_data == src.nodata, np.nan, ecostress_data)
ecostress_transform = src.transform
return ecostress_data, ecostress_transform
if __name__ == "__main__":
# Fetching Data Progress Bar
with tqdm(total=100, desc='Fetching Data', unit='%', ncols=80) as pbar:
# Fetch Traxelektronik data
trax_stations, trax_temperatures = fetch_traxelektronik_data()
pbar.update(20)
# Geocode Traxelektronik stations
geolocator = Nominatim(user_agent="geoapi")
trax_coords = {}
for station in tqdm(trax_stations, desc="Geocoding Traxelektronik Stations", unit="station"):
cleaned_station = clean_station_name(station)
coords, was_cleaned = get_coordinates(cleaned_station, geolocator)
if coords:
trax_coords[cleaned_station] = coords # Use cleaned station names as keys
pbar.update(30)
# Fetch Netatmo data
netatmo_stations, netatmo_temperatures, netatmo_coords = fetch_netatmo_data()
pbar.update(20)
# Fetch IMGW data
imgw_stations, imgw_temperatures, imgw_station_ids = fetch_imgw_data()
pbar.update(10)
# Geocode IMGW stations
imgw_coords = {}
for station in tqdm(imgw_stations, desc="Geocoding IMGW Stations", unit="station"):
coords, was_cleaned = get_coordinates(station, geolocator)
if coords:
imgw_coords[station] = coords
pbar.update(20)
# Paths to MODIS data folders
modis_day_folder = 'input/MODIS_day'
modis_night_folder = 'input/MODIS_night'
# Call the function with all required arguments
smoothed_temperature, transform = plot_stations_with_boundary_and_interpolation(
trax_stations, trax_temperatures, trax_coords,
netatmo_stations, netatmo_temperatures, netatmo_coords,
imgw_stations, imgw_temperatures, imgw_coords,
modis_day_folder, modis_night_folder
)
# Coordinates to test
test_lat = 51.40118
test_lon = 21.14805
interpolated_temp = get_interpolated_temperature_at_point(
test_lat, test_lon, smoothed_temperature, transform
)
if interpolated_temp is not None:
print(f"\nInterpolated temperature value at {test_lat}, {test_lon} is {interpolated_temp:.2f}°C")
else:
print(f"\nCoordinates {test_lat}, {test_lon} are outside the interpolated grid or have no data.")
Here, can you also additionally fix the color scale? It's definitely not that one in my input/color_scale.csv, if anything it is written in this format:
value,color
-40,#ff6eff
-38,#ff46f8
-37,#f627eb
etc..
But they're definitely not representing the current values exactly. They're less contrast, etc.ASSISTANT
Certainly! I can help you fix the color scale in your script so that it correctly reflects the values in your `input/color_scale.csv` file with proper contrast.
**Problem Analysis:**
- The issue seems to be that the color mapping in the plot is not correctly representing the values in your `color_scale.csv` file.
- This might be due to the way the colormap is being created and normalized in the code.
**Possible Causes:**
1. **Normalization Range Mismatch:**
- In your script, the normalization range (`vmin` and `vmax`) is set based on the minimum and maximum values from the `color_scale.csv` file.
- If your data (e.g., temperatures) only cover a small portion of that range (e.g., from 5°C to 25°C), but the color scale ranges from -40°C to 50°C, then most of the color map will not be utilized, resulting in less contrast.
2. **Incorrect Color Mapping:**
- The way the colors are assigned to the temperature values might not be properly aligned, causing the plot to not reflect the colors as intended.
**Solution:**
- Adjust the code to create a colormap that covers only the range of temperature values present in your data, ensuring the color scale correctly reflects the data range and improves contrast.
- Ensure that the normalization (`norm`) is set to the data range, not the full range from the color scale file.
**Code Modifications:**
1. **Compute the Data Range:**
```python
# Compute min and max temperature from data
min_temp = np.floor(np.nanmin(temperature_with_ecostress))
max_temp = np.ceil(np.nanmax(temperature_with_ecostress))
```
2. **Subset the Color Scale DataFrame:**
```python
# Subset the color scale to include only the values within data range
color_scale_df = color_scale_df[(color_scale_df['value'] >= min_temp) & (color_scale_df['value'] <= max_temp)]
color_scale_df = color_scale_df.sort_values('value')
values = color_scale_df['value'].values
hex_colors = color_scale_df['color'].values
```
3. **Normalize the Values:**
```python
# Normalize the values to 0..1 over the data range
normalized_values = (values - min_temp) / (max_temp - min_temp)
```
4. **Create the Colormap:**
```python
# Create a list of tuples (position, color)
color_tuples = list(zip(normalized_values, hex_colors))
# Create a LinearSegmentedColormap
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', color_tuples)
```
5. **Set the Normalization:**
```python
# Set the normalization to the data range
norm = mcolors.Normalize(vmin=min_temp, vmax=max_temp)
```
6. **Update the Colorbar Ticks:**
```python
# Adjust colorbar ticks to match data range
ticks = np.arange(min_temp, max_temp + 1, 5) # Adjust the step size as needed
cbar1.set_ticks(ticks)
cbar1.set_ticklabels(ticks.astype(int))
```
7. **Remove Unnecessary Color Code Adjustments:**
- Since your colors in the CSV file already start with `#`, you may not need to add `#` to them.
```python
# Comment out or remove this line if not needed
# hex_colors = [color if color.startswith('#') else f"#{color}" for color in hex_colors]
```
**Step-by-Step Implementation:**
Let's integrate these changes into your script.
---
### **Updated Code Section:**
**1. Reading and Processing the Color Scale:**
Place this code where you read and process the color scale, before plotting:
```python
# Read the custom color scale from the CSV file
color_scale_df = pd.read_csv('input/color_scale.csv')
# Compute min and max temperature from data
min_temp = np.floor(np.nanmin(temperature_with_ecostress))
max_temp = np.ceil(np.nanmax(temperature_with_ecostress))
# Subset the color scale to include only the values within data range
color_scale_df = color_scale_df[(color_scale_df['value'] >= min_temp) & (color_scale_df['value'] <= max_temp)]
color_scale_df = color_scale_df.sort_values('value')
values = color_scale_df['value'].values
hex_colors = color_scale_df['color'].values
# Normalize the values to 0..1 over the data range
normalized_values = (values - min_temp) / (max_temp - min_temp)
# Create a list of tuples (position, color)
color_tuples = list(zip(normalized_values, hex_colors))
# Create a LinearSegmentedColormap
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', color_tuples)
# Set the normalization to the data range
norm = mcolors.Normalize(vmin=min_temp, vmax=max_temp)
```
**Explanation:**
- The script reads your `color_scale.csv` file.
- It computes the minimum and maximum temperature values from your `temperature_with_ecostress` data.
- It filters the color scale to include only the values within your data range.
- It normalizes these values over the range `[0, 1]` based on your data range.
- It creates a custom colormap using the filtered and normalized color values.
- It sets the normalization to your data range for plotting.
**2. Plotting with the New Color Scale:**
In your plotting code, when you use `imshow`, ensure you pass the `cmap` and `norm` variables:
```python
# Plotting Temperature Map with ECOSTRESS Adjustments
fig1, ax1 = plt.subplots(figsize=(15, 15))
# Plot the interpolated temperature grid with ECOSTRESS adjustments
img1 = ax1.imshow(temperature_with_ecostress, extent=extent, origin='upper',
cmap=cmap, norm=norm, alpha=1)
```
**3. Adjusting Colorbar Ticks:**
After creating the colorbar, adjust the ticks to match your data range:
```python
# Adjust colorbar position and style
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="3%", pad=0.1)
cbar1 = fig1.colorbar(img1, cax=cax)
# Adjust colorbar ticks to match data range
ticks = np.arange(min_temp, max_temp + 1, 5) # Adjust the step size as needed
cbar1.set_ticks(ticks)
cbar1.set_ticklabels(ticks.astype(int))
cbar1.ax.tick_params(labelsize=14)
```
**4. Remove Unnecessary Color Code Adjustments:**
Since your color codes already start with `#`, you don't need to adjust them. Comment out or remove this line:
```python
# Ensure colors start with '#'
# hex_colors = [color if color.startswith('#') else f"#{color}" for color in hex_colors]
```
**Full Implementation Example:**
Here's how the updated code section would look in your script:
```python
# Read the custom color scale from the CSV file
color_scale_df = pd.read_csv('input/color_scale.csv')
# Compute min and max temperature from data
min_temp = np.floor(np.nanmin(temperature_with_ecostress))
max_temp = np.ceil(np.nanmax(temperature_with_ecostress))
# Subset the color scale to include only the values within data range
color_scale_df = color_scale_df[(color_scale_df['value'] >= min_temp) & (color_scale_df['value'] <= max_temp)]
color_scale_df = color_scale_df.sort_values('value')
values = color_scale_df['value'].values
hex_colors = color_scale_df['color'].values
# Normalize the values to 0..1 over the data range
normalized_values = (values - min_temp) / (max_temp - min_temp)
# Create a list of tuples (position, color)
color_tuples = list(zip(normalized_values, hex_colors))
# Create a LinearSegmentedColormap
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', color_tuples)
# Set the normalization to the data range
norm = mcolors.Normalize(vmin=min_temp, vmax=max_temp)
# Plotting Temperature Map with ECOSTRESS Adjustments
fig1, ax1 = plt.subplots(figsize=(15, 15))
# Plot the interpolated temperature grid with ECOSTRESS adjustments
img1 = ax1.imshow(temperature_with_ecostress, extent=extent, origin='upper',
cmap=cmap, norm=norm, alpha=1)
# ... (Rest of your plotting code)
# Adjust colorbar position and style
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="3%", pad=0.1)
cbar1 = fig1.colorbar(img1, cax=cax)
# Adjust colorbar ticks to match data range
ticks = np.arange(min_temp, max_temp + 1, 5) # Adjust the step size as needed
cbar1.set_ticks(ticks)
cbar1.set_ticklabels(ticks.astype(int))
cbar1.ax.tick_params(labelsize=14)
# ... (Rest of your plotting code)
```
**Additional Notes:**
- **Handling Missing Values in Color Scale:**
- Ensure that your `color_scale.csv` file includes all possible temperature values in your data range.
- If not, you might want to interpolate colors for missing values or extend the range as needed.
- **Adjusting the Colorbar Step Size:**
- The `np.arange` function in the ticks setting allows you to specify the step size.
- Adjust `5` to a smaller number (e.g., `2` or `1`) if the temperature range is small.
- **Maintain Contrast:**
- By aligning the colormap and normalization range with your data range, the colors in your plot should now accurately represent the temperature variations, enhancing the contrast.
**Conclusion:**
These changes should ensure that your color scale accurately reflects the data you're plotting, improving the contrast and making the colors represent the temperature values as intended.
**Example Output:**
With these changes, the temperature map you generate should now have a color gradient that accurately reflects the temperature data, using the colors specified in your `color_scale.csv` file.
**Let me know if you have any further questions or need assistance with implementing these changes.**