import gemgis as gg
import geopandas as gpd

# Load vector data
gdf = gpd.read_file('geology.shp')

# Extract XYZ from geometry with elevation from DEM
interfaces = gg.vector.extract_xyz(gdf=gdf, dem='dem.tif')

# Define model extent
extent = [x_min, x_max, y_min, y_max]

# Clip to extent
gdf_clipped = gg.vector.clip_by_extent(gdf=gdf, extent=extent)

contacts = gpd.read_file('contacts.shp')
interfaces = gg.vector.extract_xyz(gdf=contacts, dem='dem.tif')
interfaces['formation'] = contacts['formation']
# Returns DataFrame with X, Y, Z, formation columns

measurements = gpd.read_file('structural_measurements.shp')
orientations = gg.vector.extract_xyz(gdf=measurements, dem='dem.tif')
orientations['dip'] = measurements['dip']
orientations['azimuth'] = measurements['strike'] + 90  # strike to dip direction
orientations['formation'] = measurements['formation']

profile = gg.raster.sample_from_raster(
    raster='dem.tif',
    line=[(500000, 5600000), (510000, 5605000)],
    n_samples=100
)
# Returns dict with 'distance' and 'Z' arrays

# From coordinates
extent = gg.utils.set_extent(
    x_min=500000, x_max=510000,
    y_min=5600000, y_max=5610000
)

# From GeoDataFrame bounds
extent = gg.utils.set_extent_from_bounds(gdf)

from shapely.geometry import box
extent_poly = box(500000, 5600000, 510000, 5610000)
gdf_clipped = gdf.clip(extent_poly)

# Or using gemgis
gdf_clipped = gg.vector.clip_by_extent(
    gdf=gdf,
    extent=[500000, 510000, 5600000, 5610000]
)

# Always work in projected CRS (meters) for modeling
gdf_utm = gdf.to_crs('EPSG:32632')  # UTM zone 32N

# Or use GemGIS utility
gdf_utm = gg.vector.reproject(gdf, 'EPSG:32632')