Turning KML into point shapefile
import os
import zipfile
from pathlib import Path
import geopandas as gpd
import rasterio
from rasterio.features import rasterize
from shapely.geometry import box
import numpy as np
# Paths
kmz_dir = Path("D:/RG_Inventory/Google_earth") # Directory containing KMZ files
output_shapefiles_dir = Path("D:/RG_Inventory/shapefile") # Output for shapefiles
composite_raster_path = Path("D:/RG_Inventory/composite/composite_raster.tif") # Final raster
boundary_raster_path = "D:/october_350/masked_3500.tif" # Boundary raster
output_shapefiles_dir.mkdir(parents=True, exist_ok=True)
# Read the boundary raster to define the extent and resolution
with rasterio.open(boundary_raster_path) as src:
boundary_transform = src.transform
boundary_crs = src.crs
boundary_shape = (src.height, src.width) # (rows, columns)
boundary_bounds = src.bounds
# Initialize composite data array with zeros
composite_data = np.zeros(boundary_shape, dtype="uint8")
# Process each KMZ file
for kmz_file in kmz_dir.glob("*.kmz"):
try:
# Extract KMZ
extracted_dir = output_shapefiles_dir / kmz_file.stem
extracted_dir.mkdir(exist_ok=True)
with zipfile.ZipFile(kmz_file, 'r') as z:
z.extractall(extracted_dir)
# Locate KML file
kml_files = list(extracted_dir.rglob("*.kml"))
if not kml_files:
print(f"No KML file found in {kmz_file}")
continue
kml_file = kml_files[0] # Assuming one KML per KMZ
# Read KML as GeoDataFrame
point_gdf = gpd.read_file(kml_file, driver='KML')
# Check CRS of the points and reproject to match the raster CRS
if point_gdf.crs != boundary_crs:
print(f"Reprojecting {kmz_file.name} from {point_gdf.crs} to {boundary_crs}")
point_gdf = point_gdf.to_crs(boundary_crs)
point_gdf['label'] = 1 # Mark points as 1
# Save as individual shapefile
output_shapefile = output_shapefiles_dir / f"{kmz_file.stem}.shp"
point_gdf.to_file(output_shapefile, driver="ESRI Shapefile")
print(f"Saved shapefile: {output_shapefile}")
# Check if points are within the boundary bounds
points_within_bounds = point_gdf.cx[boundary_bounds[0]:boundary_bounds[2],
boundary_bounds[1]:boundary_bounds[3]]
if points_within_bounds.empty:
print(f"No points within boundary for {kmz_file.name}")
# Rasterize points into the composite raster
shapes = ((geom, 1) for geom in point_gdf.geometry)
rasterized = rasterize(
shapes=shapes,
out_shape=boundary_shape,
transform=boundary_transform,
fill=0,
dtype="uint8"
composite_data += rasterized # Add to the composite raster
except Exception as e:
print(f"Error processing {kmz_file.name}: {e}")
# Save the composite raster
with rasterio.open(
composite_raster_path,
"w",
driver="GTiff",
height=boundary_shape[0],
width=boundary_shape[1],
count=1,
dtype="uint8",
crs=boundary_crs,
transform=boundary_transform,
) as dst:
dst.write(composite_data, 1)
print(f"Composite raster saved to {composite_raster_path}")
Overlay points giving them 0, 1
import geopandas as gpd
import rasterio
import numpy as np
from rasterio.features import geometry_mask
# Step 1: Open the raster (ROI) file to get its metadata and the extent
raster_path = "D:/october_350/masked_3500.tif"
with rasterio.open(raster_path) as src:
# Get the raster's metadata
meta = src.meta
# Get the raster's affine transformation
transform = src.transform
# Read the raster data (optional, in case you want to visualize it)
raster_data = src.read(1) # assuming it's a single-band raster
# Step 2: Prepare an empty array of the same shape as the raster filled with 0
output_array = np.zeros_like(raster_data)
# Step 3: Loop through each shapefile and update the output_array
for shapefile_path in ["D:/RG_Inventory/shapefile/RG 12.shp", "D:/RG_Inventory/shapefile/RG.shp",
"D:/RG_Inventory/shapefile/rg3.shp", "D:/RG_Inventory/shapefile/RG4.shp",
"D:/RG_Inventory/shapefile/RG5.shp", "D:/RG_Inventory/shapefile/rg6.shp",
"D:/RG_Inventory/shapefile/RG7.shp", "D:/RG_Inventory/shapefile/RG8.shp",
"D:/RG_Inventory/shapefile/RG9.shp", "D:/RG_Inventory/shapefile/RG10.shp"]:
# Read the shapefile
shapefile = gpd.read_file(shapefile_path)
# Create a mask for the current shapefile
mask = geometry_mask(shapefile.geometry, transform=transform, invert=True,
out_shape=(meta['height'], meta['width']))
# Add the shapefile's values to the output array (set corresponding mask areas to 1)
output_array[mask] = 1
# Step 4: Save the resulting raster
output_raster_path = "D:/RG_Inventory/composite/output_raster.tif"
with rasterio.open(output_raster_path, 'w', **meta) as dst:
dst.write(output_array, 1) # Write the output array to the first band
