EDeMOS/basegridhex.py at main · EnergySystemsModellingLab/EDeMOS · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
"""BaseGridHex.ipynb


# Part 1. Create base grid with H3

### Import necessary modules
"""

# Spatial
import geopandas as gpd
from geopandas.tools import sjoin

# Mapping / Plotting
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar

import config

"""### Functions for creating heaxgons"""

from create_hex import*

"""### Define area of interest"""

area = config.AREA_OF_INTEREST
print(area)


"""### Import layers to be used"""

## admininstrative boundary
if area == "COUNTRY":
    admin_gdf = gpd.read_file(config.ADMIN_PATH / config.ADMIN_GPKG, layer=config.ADMIN_LAYER_COUNTRY)
    region_gdf = gpd.read_file(config.ADMIN_PATH / config.ADMIN_GPKG, layer=config.ADMIN_LAYER_REGION)
else:
    region_gdf = gpd.read_file(config.ADMIN_PATH / config.ADMIN_GPKG, layer=config.ADMIN_LAYER_REGION)
    region_gdf = region_gdf[region_gdf[config.ADMIN_REGION_COLUMN_NAME]==area]
    admin_gdf = region_gdf

print(admin_gdf.crs)

"""### H3 - Hexagon - grid"""

print("Creating a buffer to ensure full hexagon coverage...")
# Define buffer distance in meters.
buffer_distance_meters = config.buffer_distance_meters

# Store original CRS
original_crs = admin_gdf.crs
# Reproject to a projected CRS (e.g., a UTM zone) for accurate buffering in meters.
admin_gdf_proj = admin_gdf.to_crs(config.CRS_PROJ)

# Create a single unified geometry for buffering.
unified_geometry = admin_gdf_proj.union_all()
# Apply the buffer
buffered_geometry_proj = unified_geometry.buffer(buffer_distance_meters)

# Create a new GeoDataFrame for the buffered area
admin_gdf_buffered_proj = gpd.GeoDataFrame(geometry=[buffered_geometry_proj], crs=config.CRS_PROJ)
# Reproject the buffered GeoDataFrame back to the original CRS (WGS84)
admin_gdf_buffered = admin_gdf_buffered_proj.to_crs(original_crs)

print("Buffer created successfully.")

size = config.HEX_SIZE ## resolution info here https://h3geo.org/docs/core-library/restable
# hexagons = feat(admin_gdf, size)
hexagons_unclipped = feat(admin_gdf_buffered, size)
print("Clipping hexagons and attaching region attributes...")
hexagons = gpd.sjoin(hexagons_unclipped,  region_gdf[[config.ADMIN_REGION_COLUMN_NAME, "geometry"]], how="inner", predicate="intersects")
# The sjoin adds an 'index_right' column
hexagons = hexagons.drop(columns=['index_right'])
hexagons = hexagons.drop(columns=['index'])
hexagons = hexagons.drop_duplicates(subset='h3_index').reset_index(drop=True)
hexagons['id'] = range(1, len(hexagons)+1)
print(hexagons.columns)
# is_unique = hexagons['h3_index'].is_unique
# print(f"Are all h3_index values unique? {is_unique}")
# hexagons.set_index('h3_index')

"""#### Select base map grid"""
def plot_and_save_map(hexagons, admin_gdf, region_gdf):
    # Plot basemap
    plt.rcParams.update({'font.size': 22})
    fig, ax = plt.subplots(figsize=(25, 15))

    # Convert the dataset to a coordinate
    hexagons.plot(ax=ax, edgecolor='brown', alpha=0.2)
    admin_gdf.plot(ax=ax, edgecolor='brown', alpha=0.2)
    region_gdf.plot(ax=ax, edgecolor='brown', alpha=0.2)
    ax.set_aspect('equal', 'box')
    # Add latitude and longitude labels
    ax.set_xlabel('Longitude (°)')
    ax.set_ylabel('Latitude (°)')

    # Compute the distance-per-pixel of the map
    # see https://geopandas.org/en/latest/gallery/matplotlib_scalebar.html#Geographic-coordinate-system-(degrees)
    assert admin_gdf.crs == config.CRS_WGS84
    from shapely.geometry.point import Point
    points = gpd.GeoSeries(
        [Point(-73.5, 40.5), Point(-74.5, 40.5)], crs=4326
    )  # Geographic WGS 84 - degrees
    points = points.to_crs(32619)  # Projected WGS 84 - meters
    distance_meters = points[0].distance(points[1])

    # Add a scale bar
    scalebar = ScaleBar(
        distance_meters,
        dimension="si-length",
        location='lower left',
        length_fraction=0.1,
        width_fraction=0.001,
        units='m',
        color='black',
        fixed_value=None
    )

    ax.add_artist(scalebar)

    # plt.show()
    # Save plot as figure
    plt.savefig(config.OUTPUT_DIR / f'admin_level_basemap_{config.COUNTRY}.png', bbox_inches='tight')
    plt.savefig(config.OUTPUT_DIR / f'admin_level_basemap_{config.COUNTRY}.svg', format='svg', bbox_inches='tight')
    print(f"Map saved to {config.OUTPUT_DIR}")

plot_and_save_map(hexagons, admin_gdf, region_gdf)

# Export dataframe to csv or gpkg
#hexagons.to_csv(config.OUTPUT_DIR  + "\\" + f'h3_grid_at_hex_{size}.csv', index=False)
# hexagons.to_file(config.OUTPUT_DIR / f'h3_grid_at_hex_{size}.shp', index=False)
hexagons.to_file(config.OUTPUT_DIR / config.H3_GRID_HEX_SHP, index=False) # file used in the other scripts
# hexagons.to_file(config.OUTPUT_DIR / "hex.geojson")
# admin_gdf.to_file(config.OUTPUT_DIR / f'area_gdf.gpkg', index=False)
# admin_gdf.to_file(config.OUTPUT_DIR  / f'area_gdf.geojson', driver='GeoJSON', index=False)