-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathvisualize_cumulative_points_with_counts.py
55 lines (43 loc) · 1.98 KB
/
visualize_cumulative_points_with_counts.py
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
import os
import geopandas as gpd
from datetime import datetime, timedelta
import json
from shapely.geometry import shape
# Define directories
base_dir = '.'
cumulative_dir = os.path.join(base_dir, 'cumulative')
shapefile_path = os.path.join(base_dir, 'basisdaten', 'VG5000_GEM.shp')
output_dir = os.path.join(base_dir, 'shapefile_cumulative')
# Create output directory if it doesn't exist
os.makedirs(output_dir, exist_ok=True)
# Load the shapefile
shapefile_gdf = gpd.read_file(shapefile_path)
shapefile_gdf.set_crs(epsg=25832, inplace=True) # Set initial CRS to EPSG:25832
shapefile_gdf = shapefile_gdf.to_crs(epsg=3857) # Transform to match the map's CRS
# Introduce the overwrite variable
overwrite = False
# Process each cumulative file
for cumulative_file in sorted(os.listdir(cumulative_dir)):
if cumulative_file.endswith('.geojson'):
cumulative_path = os.path.join(cumulative_dir, cumulative_file)
date = cumulative_file.split('_')[0]
output_file_path = os.path.join(output_dir, f'{date}_VG5000_GEM_with_counts.shp')
# Check if the output file exists and overwrite is False
if os.path.exists(output_file_path) and not overwrite:
continue
# Load the GeoJSON file
with open(cumulative_path, 'r') as f:
geojson_data = json.load(f)
# Create a GeoDataFrame for the GeoJSON points
points_gdf = gpd.GeoDataFrame.from_features(geojson_data['features'])
points_gdf.set_crs(epsg=4326, inplace=True)
points_gdf = points_gdf.to_crs(epsg=3857)
# Initialize the 'NUMPOINTS' column
shapefile_gdf['NUMPOINTS'] = 0
# Count points within each polygon
for idx, poly in shapefile_gdf.iterrows():
count = points_gdf.within(poly.geometry).sum()
shapefile_gdf.at[idx, 'NUMPOINTS'] = count
# Save the updated shapefile
shapefile_gdf.to_file(output_file_path)
print(f'Processed {cumulative_file} and saved to {output_file_path}')