V. Raster Analysis

import geopandas as gpd
import pandas as pd
import osmnx as ox
import matplotlib.pyplot as plt
import numpy as np
import rasterio

#Можно выгрузить данные с портала World Pop (300 Mb + )

# # Ссылка на данные WorldPop
# url = "https://data.worldpop.org/GIS/Population/Global_2000_2020_Constrained/2020/BSGM/RUS/rus_ppp_2020_constrained.tif"

# # Имя файла для сохранения
# output_file = "worldpop_russia_2020.tif"

# # Скачивание файла
# response = requests.get(url, stream=True)

# if response.status_code == 200:
#     with open(output_file, "wb") as file:
#         for chunk in response.iter_content(chunk_size=1024):
#             file.write(chunk)
#     print(f"File saved as {output_file}")
# else:
#     print("Failed to download file. Status code:", response.status_code)

output_file = './data/tula_region_population.tif'


with rasterio.open(output_file) as dataset:
    print("CRS:", dataset.crs)
    print("Bounds:", dataset.bounds)
    print("Resolution:", dataset.res)

with rasterio.open(output_file) as src:
    data = src.read(1)  

    plt.figure(figsize=(10, 10))
    plt.imshow(data, cmap='viridis') 
    plt.colorbar(label="Population Density") 
    plt.title("Population Density Tula Region (WorldPop, Russia, 2020)")
    plt.show()

with rasterio.open(output_file) as src:
    data = src.read(1) 
    # Выводим статистику
    print("Min value:", np.min(data))
    print("Max value:", np.max(data))
    print("Mean value:", np.mean(data))

with rasterio.open(output_file) as src:
    data = src.read(1, masked=True)  # NoData будет маской

print("Min:", data.min())
print("Max:", data.max())
print("Mean:", data.mean())

with rasterio.open(output_file) as src:
    data = src.read(1, masked=True)

plt.figure(figsize=(10, 10))
plt.imshow(data, cmap="viridis")
plt.colorbar(label="Population Density")
plt.title("Population Density Tula Region (WorldPop, Russia, 2020)")
plt.show()

area_name = "Тула, Россия"

gdf_borders = ox.geocode_to_gdf(area_name)

gdf_borders.explore(tiles='cartodbpositron')

from rasterio.mask import mask
import json

# Убедимся, что CRS совпадает
with rasterio.open(output_file) as src:
    
    # Если CRS различаются — перепроецируем границы
    if gdf_borders.crs != src.crs:
        gdf_borders = gdf_borders.to_crs(src.crs)
    
    # Преобразуем геометрию в формат GeoJSON
    geometries = [json.loads(gdf_borders.to_json())['features'][0]['geometry']]
    
    # Обрезка по точной геометрии
    out_image, out_transform = mask(src, geometries, crop=True)
    
    # Обновляем профиль
    out_meta = src.meta.copy()
    out_meta.update({
        "driver": "GTiff",
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    })

file_cropped_path = "./data/tula_cropped_population.tif"

# Сохраняем результат
with rasterio.open(file_cropped_path, "w", **out_meta) as dest:
    dest.write(out_image)

with rasterio.open(file_cropped_path) as src:
    # Вывод основной информации о растре
    print("CRS:", src.crs)  # Система координат
    print("Bounds:", src.bounds)  # Границы растра
    print("Width, Height:", src.width, src.height)  # Размер растра
    print("Number of bands:", src.count)  # Количество слоев
    print("Data type:", src.dtypes)  # Тип данных
    print("Transform:", src.transform)  # Аффинная трансформация

with rasterio.open(file_cropped_path) as src:
    transform = src.transform  # Аффинная трансформация растра

    # Извлечение размеров пикселя
    pixel_width = transform.a  # Размер пикселя по оси X
    pixel_height = -transform.e  # Размер пикселя по оси Y (берем с минусом, так как Y направлен вниз в системе координат)

    print(f"Pixel Width: {pixel_width} units")
    print(f"Pixel Height: {pixel_height} units")

target_crs = gdf_borders.estimate_utm_crs()

from rasterio.warp import calculate_default_transform, reproject, Resampling

output_file_cropped = "./data/tula_cropped_reprojected.tif"

with rasterio.open(file_cropped_path) as src:
    transform, width, height = calculate_default_transform(src.crs, target_crs, src.width, src.height, *src.bounds)
    profile = src.profile
    profile.update({
        'crs': target_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rasterio.open(output_file_cropped, 'w', **profile) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=target_crs,
                resampling=Resampling.nearest
            )

with rasterio.open(output_file_cropped) as src:
    transform = src.transform 

    # Извлечение размеров пикселя
    pixel_width = transform.a  # Размер пикселя по оси X
    pixel_height = -transform.e  # Размер пикселя по оси Y (берем с минусом, так как Y направлен вниз в системе координат)

    print(f"Pixel Width: {pixel_width} units")
    print(f"Pixel Height: {pixel_height} units")

with rasterio.open(output_file_cropped) as src:
    data = src.read(1, masked=True)  

    plt.figure(figsize=(10, 10))
    plt.imshow(data, cmap='viridis') 
    plt.colorbar(label="Population Density")  
    plt.title("Population Density (WorldPop, Russia, 2020)")
    plt.show()

from shapely.geometry import box

def create_regular_grid(gdf, square_size):
    # Если набор в географических координатах, перепроецируем в UTM
    utm_crs = gdf.estimate_utm_crs()
    gdf_utm = gdf.to_crs(utm_crs)
    # Получаем границы (bounding box) области
    minX, minY, maxX, maxY = gdf_utm.total_bounds
    grid_cells = []
    x, y = minX, minY
    # Двойной цикл: заполняем прямоугольник границ квадратами заданного размера
    while y <= maxY:
        while x <= maxX:
            grid_cells.append(box(x, y, x + square_size, y + square_size))
            x += square_size # двигаемся по X
        x = minX             # возврат к началу по X
        y += square_size     # двигаемся по Y (следующий ряд)
    # Формируем GeoDataFrame
    fishnet = gpd.GeoDataFrame(geometry=grid_cells, crs=utm_crs)
    return fishnet


# Строим сетку с шагом 1000 м
fishnet = create_regular_grid(gdf_borders, square_size=1000)
print(f"Создано ячеек: {len(fishnet)}")


fishnet.explore(tiles='cartodbpositron')

print(fishnet.crs)
print(target_crs)

import rasterstats as rs

# Рассчитываем зональную статистику (сумма значений растра для каждого полигона)
stats = rs.zonal_stats(fishnet, output_file_cropped, stats="sum", geojson_out=True)

# Преобразуем результат в GeoDataFrame
gdf_stats = gpd.GeoDataFrame.from_features(stats, crs=target_crs)

# Теперь у вас есть GeoDataFrame с дополнительной колонкой, содержащей сумму значений для каждого полигона
gdf_stats.head()

gdf_stats['density'] = gdf_stats['sum']/(gdf_stats.geometry.area/1000000)

gdf_stats.explore(column='density', tiles='cartodbpositron')

mkd_df = pd.read_csv('data/tula_mkd_geocoded.csv')

mkd_gdf = gpd.GeoDataFrame(mkd_df, geometry=gpd.points_from_xy(mkd_df.lon, mkd_df.lat), crs="EPSG:4326")

#mkd_gdf.explore(tiles='cartodbpositron')

# Приводим жилую площадь к числовому формату:
# заменяем запятую на точку и конвертируем в float
mkd_gdf["area_residential"] = (
    mkd_gdf["area_residential"]
        .str.replace(",", ".", regex=False)
)

mkd_gdf["area_residential"] = pd.to_numeric(
    mkd_gdf["area_residential"],
    errors="coerce"   # некорректные значения -> NaN
)

# Оценка численности населения (28 м² на человека - приблизительное число)
mkd_gdf['estimated_population'] = mkd_gdf['area_residential'] / 28

# Смотрим на результат
mkd_gdf.head()

# Приводим слой МКД к целевой системе координат (CRS сетки)
mkd_gdf_utm = mkd_gdf.to_crs(target_crs)

# Пространственное соединение:
# определяем, в какую ячейку сетки попадает каждый дом
msk_in_grid = gpd.sjoin(mkd_gdf_utm, gdf_stats, predicate='within')

# Агрегирование: суммируем оценочную численность населения по ячейкам сетки
# 'index_right' — индекс ячейки, в которую попала точка
pop_grid_df = (
    msk_in_grid
    .groupby("index_right")["estimated_population"]
    .sum()
    .reset_index(name="pop_sum")
)

# Присоединяем агрегированные значения к слою сетки
pop_grid_gdf = gdf_stats.join(
    pop_grid_df.set_index("index_right")["pop_sum"]
)

# Если в ячейке отсутствуют МКД, заменяем NaN на 0
pop_grid_gdf["pop_sum"] = pop_grid_gdf["pop_sum"].fillna(0)

# Расчёт плотности населения (чел/км²):
# делим численность на площадь ячейки (в км²)
pop_grid_gdf["pop_density_mkd"] = (
    pop_grid_gdf["pop_sum"] / (pop_grid_gdf.geometry.area / 1_000_000)
)

pop_grid_gdf.explore(column='pop_density_mkd', cmap='YlGnBu', tiles='cartodbpositron', scheme='NaturalBreaks', k=5, legend=True, missing_kwds={'color': '#ffffff00','fillOpacity': 0})

pop_grid_gdf["density_diff"] = pop_grid_gdf["pop_density_mkd"] - pop_grid_gdf["density"]

pop_grid_gdf["density_diff"].describe()

max_abs = max(
    abs(pop_grid_gdf["density_diff"].min()),
    abs(pop_grid_gdf["density_diff"].max())
)

pop_grid_gdf.explore(
    tiles='cartodbpositron',
    column="density_diff",
    cmap="coolwarm",
    vmin=-max_abs,
    vmax= max_abs,
    tooltip=["density_diff", "pop_density_mkd", "density"],
    legend=True,
    missing_kwds={'color': '#ffffff00','fillOpacity': 0}
)