package service import ( "bugulma/backend/internal/domain" "bugulma/backend/internal/geospatial" "context" "fmt" "math" "gorm.io/gorm" ) // GeospatialService provides advanced geospatial operations for sites and geographical features type GeospatialService struct { db *gorm.DB geoFeatureRepo domain.GeographicalFeatureRepository } // NewGeospatialService creates a new geospatial service func NewGeospatialService(db *gorm.DB, geoFeatureRepo domain.GeographicalFeatureRepository) *GeospatialService { return &GeospatialService{ db: db, geoFeatureRepo: geoFeatureRepo, } } // SpatialQuery represents a spatial query with various criteria type SpatialQuery struct { CenterLat float64 CenterLng float64 RadiusKm float64 SiteTypes []domain.SiteType ResourceTypes []domain.ResourceType MaxResults int } // SpatialResult represents a site with spatial metadata type SpatialResult struct { Site *domain.Site DistanceKm float64 Bearing float64 // Direction from center point DrivingTime *int // Estimated driving time in minutes (optional) } // SpatialCluster represents a cluster of spatially close sites type SpatialCluster struct { ID int `json:"id"` CentroidLat float64 `json:"centroid_lat"` CentroidLng float64 `json:"centroid_lng"` SiteCount int `json:"site_count"` Sites []SpatialResult `json:"sites"` RadiusKm float64 `json:"radius_km"` } // FindNearbySites finds sites within radius with advanced filtering func (gs *GeospatialService) FindNearbySites(ctx context.Context, query SpatialQuery) ([]SpatialResult, error) { var results []SpatialResult // Build the base query with PostGIS using GeoHelper to centralize PostGIS fragments geo := geospatial.NewGeoHelper(gs.db) baseQuery := ` SELECT s.*, ST_Distance(s.location_geometry::geography, ` + geo.PointExpr() + `::geography) / 1000 as distance_km, ST_Azimuth(` + geo.PointExpr() + `, s.location_geometry) as bearing FROM sites s WHERE s.location_geometry IS NOT NULL AND ` + geo.DWithinExpr("s.location_geometry") + ` ` // Placeholders appear in this order in the query: distance (lng,lat), azimuth (lng,lat), dwithin point (lng,lat), radius args := []interface{}{query.CenterLng, query.CenterLat, query.CenterLng, query.CenterLat, query.CenterLng, query.CenterLat, query.RadiusKm} // Add site type filter if len(query.SiteTypes) > 0 { placeholders := "" for i, siteType := range query.SiteTypes { if i > 0 { placeholders += "," } placeholders += "?" args = append(args, siteType) } baseQuery += fmt.Sprintf(" AND s.site_type IN (%s)", placeholders) } // Order by distance and limit results // Add ORDER BY distance using the parameterized point again baseQuery += " ORDER BY s.location_geometry <-> " + geo.PointExpr() // For ORDER BY we need to append the final lng/lat placeholders args = append(args, query.CenterLng, query.CenterLat) if query.MaxResults > 0 { baseQuery += " LIMIT ?" args = append(args, query.MaxResults) } // Execute query rows, err := gs.db.Raw(baseQuery, args...).Rows() if err != nil { return nil, err } defer rows.Close() for rows.Next() { var site domain.Site var distanceKm, bearing float64 var locationGeometry interface{} // Geometry column - we don't use it in the struct // Scan the row - need to be careful with column order // Note: location_geometry is scanned but not stored in the struct err := rows.Scan( &site.ID, &site.Name, &site.Latitude, &site.Longitude, &locationGeometry, &site.SiteType, &site.FloorAreaM2, &site.Ownership, &site.OwnerOrganizationID, &site.AvailableUtilities, &site.ParkingSpaces, &site.LoadingDocks, &site.CraneCapacityTonnes, &site.EnergyRating, &site.WasteManagement, &site.EnvironmentalImpact, &site.YearBuilt, &site.BuilderOwner, &site.Architect, &site.OriginalPurpose, &site.CurrentUse, &site.Style, &site.Materials, &site.Storeys, &site.HeritageStatus, &site.Notes, &site.Sources, &site.CreatedAt, &site.UpdatedAt, &distanceKm, &bearing, ) if err != nil { return nil, err } results = append(results, SpatialResult{ Site: &site, DistanceKm: distanceKm, Bearing: bearing, }) } return results, nil } // CalculateDistanceMatrix calculates distances between multiple points efficiently func (gs *GeospatialService) CalculateDistanceMatrix(ctx context.Context, points []domain.Site) ([][]float64, error) { if len(points) == 0 { return [][]float64{}, nil } matrix := make([][]float64, len(points)) for i := range matrix { matrix[i] = make([]float64, len(points)) } // For small matrices, calculate directly // For larger matrices, could use PostGIS functions for i := 0; i < len(points); i++ { for j := 0; j < len(points); j++ { if i == j { matrix[i][j] = 0 } else { matrix[i][j] = gs.calculateHaversineDistance( points[i].Latitude, points[i].Longitude, points[j].Latitude, points[j].Longitude, ) } } } return matrix, nil } // calculateHaversineDistance calculates distance using Haversine formula func (gs *GeospatialService) calculateHaversineDistance(lat1, lon1, lat2, lon2 float64) float64 { const R = 6371 // Earth radius in km dLat := (lat2 - lat1) * math.Pi / 180 dLon := (lon2 - lon1) * math.Pi / 180 a := math.Sin(dLat/2)*math.Sin(dLat/2) + math.Cos(lat1*math.Pi/180)*math.Cos(lat2*math.Pi/180)* math.Sin(dLon/2)*math.Sin(dLon/2) c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a)) return R * c } // ValidateGeometry validates PostGIS geometry func (gs *GeospatialService) ValidateGeometry(ctx context.Context, siteID string) error { var count int64 err := gs.db.Raw(` SELECT COUNT(*) FROM sites WHERE id = ? AND ST_IsValid(location_geometry) `, siteID).Scan(&count).Error if err != nil { return err } if count == 0 { return fmt.Errorf("invalid geometry for site %s", siteID) } return nil } // GetSpatialStatistics returns spatial statistics for sites func (gs *GeospatialService) GetSpatialStatistics(ctx context.Context) (map[string]interface{}, error) { var stats struct { TotalSites int64 SitesWithGeometry int64 AvgDistance float64 MaxDistance float64 MedianLatitude float64 MedianLongitude float64 } // Get basic counts gs.db.Raw("SELECT COUNT(*) as total_sites FROM sites").Scan(&stats.TotalSites) gs.db.Raw("SELECT COUNT(*) as sites_with_geometry FROM sites WHERE location_geometry IS NOT NULL").Scan(&stats.SitesWithGeometry) // Calculate spatial statistics if we have geometry data if stats.SitesWithGeometry > 1 { row := gs.db.Raw(` SELECT AVG(ST_Distance(a.location_geometry::geography, b.location_geometry::geography)) / 1000 as avg_distance, MAX(ST_Distance(a.location_geometry::geography, b.location_geometry::geography)) / 1000 as max_distance FROM sites a CROSS JOIN sites b WHERE a.id < b.id AND a.location_geometry IS NOT NULL AND b.location_geometry IS NOT NULL `).Row() row.Scan(&stats.AvgDistance, &stats.MaxDistance) // Calculate median coordinates row2 := gs.db.Raw(` SELECT PERCENTILE_CONT(0.5) WITHIN GROUP (ORDER BY latitude) as median_latitude, PERCENTILE_CONT(0.5) WITHIN GROUP (ORDER BY longitude) as median_longitude FROM sites WHERE location_geometry IS NOT NULL `).Row() row2.Scan(&stats.MedianLatitude, &stats.MedianLongitude) } return map[string]interface{}{ "total_sites": stats.TotalSites, "sites_with_geometry": stats.SitesWithGeometry, "avg_distance_km": stats.AvgDistance, "max_distance_km": stats.MaxDistance, "median_latitude": stats.MedianLatitude, "median_longitude": stats.MedianLongitude, }, nil } // FindSpatialClusters finds clusters of spatially close sites using DBSCAN algorithm func (gs *GeospatialService) FindSpatialClusters(ctx context.Context, minPoints int, radiusKm float64) ([]SpatialCluster, error) { var clusters []SpatialCluster // Use PostGIS ST_ClusterDBSCAN to group nearby sites // Convert radius from km to degrees (approximate) radiusDegrees := radiusKm / 111.32 // Rough conversion: 1 degree ≈ 111.32 km query := ` WITH clustered_sites AS ( SELECT s.*, ST_ClusterDBSCAN(s.location_geometry, ?, ?) OVER () as cluster_id, ST_Centroid(ST_Collect(s.location_geometry)) OVER (PARTITION BY ST_ClusterDBSCAN(s.location_geometry, ?, ?) OVER ()) as cluster_centroid FROM sites s WHERE s.location_geometry IS NOT NULL ), cluster_stats AS ( SELECT cluster_id, COUNT(*) as site_count, ST_Y(cluster_centroid) as centroid_lat, ST_X(cluster_centroid) as centroid_lng, MAX(ST_Distance(cluster_centroid::geography, location_geometry::geography)) / 1000 as cluster_radius_km FROM clustered_sites WHERE cluster_id IS NOT NULL GROUP BY cluster_id, cluster_centroid HAVING COUNT(*) >= ? ) SELECT cs.cluster_id, cstats.centroid_lat, cstats.centroid_lng, cstats.site_count, cstats.cluster_radius_km, cs.id as site_id, cs.name, cs.latitude, cs.longitude, cs.location_geometry, cs.site_type, cs.floor_area_m2, cs.ownership, cs.owner_organization_id, cs.available_utilities, cs.parking_spaces, cs.loading_docks, cs.crane_capacity_tonnes, cs.energy_rating, cs.waste_management, cs.environmental_impact, cs.year_built, cs.builder_owner, cs.architect, cs.original_purpose, cs.current_use, cs.style, cs.materials, cs.storeys, cs.heritage_status, cs.notes, cs.sources, cs.created_at, cs.updated_at, ST_Distance(cstats.cluster_centroid::geography, cs.location_geometry::geography) / 1000 as distance_from_centroid FROM clustered_sites cs JOIN cluster_stats cstats ON cs.cluster_id = cstats.cluster_id ORDER BY cs.cluster_id, distance_from_centroid ` rows, err := gs.db.Raw(query, radiusDegrees, minPoints, radiusDegrees, minPoints, minPoints).Rows() if err != nil { return nil, fmt.Errorf("failed to execute clustering query: %w", err) } defer rows.Close() clusterMap := make(map[int]*SpatialCluster) for rows.Next() { var ( clusterID int centroidLat float64 centroidLng float64 siteCount int clusterRadiusKm float64 site domain.Site distanceFromCentroid float64 locationGeometry interface{} // Geometry column - we don't use it in the struct ) err := rows.Scan( &clusterID, ¢roidLat, ¢roidLng, &siteCount, &clusterRadiusKm, &site.ID, &site.Name, &site.Latitude, &site.Longitude, &locationGeometry, &site.SiteType, &site.FloorAreaM2, &site.Ownership, &site.OwnerOrganizationID, &site.AvailableUtilities, &site.ParkingSpaces, &site.LoadingDocks, &site.CraneCapacityTonnes, &site.EnergyRating, &site.WasteManagement, &site.EnvironmentalImpact, &site.YearBuilt, &site.BuilderOwner, &site.Architect, &site.OriginalPurpose, &site.CurrentUse, &site.Style, &site.Materials, &site.Storeys, &site.HeritageStatus, &site.Notes, &site.Sources, &site.CreatedAt, &site.UpdatedAt, &distanceFromCentroid, ) if err != nil { return nil, fmt.Errorf("failed to scan cluster row: %w", err) } // Get or create cluster cluster, exists := clusterMap[clusterID] if !exists { cluster = &SpatialCluster{ ID: clusterID, CentroidLat: centroidLat, CentroidLng: centroidLng, SiteCount: siteCount, Sites: []SpatialResult{}, RadiusKm: clusterRadiusKm, } clusterMap[clusterID] = cluster } // Add site to cluster cluster.Sites = append(cluster.Sites, SpatialResult{ Site: &site, DistanceKm: distanceFromCentroid, Bearing: 0, // Could calculate bearing from centroid if needed }) } // Convert map to slice for _, cluster := range clusterMap { clusters = append(clusters, *cluster) } return clusters, nil } // Geographical Feature Methods // FindNearbyGeographicalFeatures finds geographical features within radius func (gs *GeospatialService) FindNearbyGeographicalFeatures(ctx context.Context, featureType domain.GeographicalFeatureType, lat, lng, radiusKm float64) ([]*domain.GeographicalFeature, error) { switch featureType { case domain.GeographicalFeatureTypeRoad: return gs.geoFeatureRepo.GetRoadsWithinRadius(ctx, lat, lng, radiusKm) case domain.GeographicalFeatureTypeGreenSpace: return gs.geoFeatureRepo.GetGreenSpacesWithinRadius(ctx, lat, lng, radiusKm) default: // For other feature types, use general spatial query return gs.findGeographicalFeaturesWithinRadius(ctx, featureType, lat, lng, radiusKm) } } // findGeographicalFeaturesWithinRadius is a helper for general feature type queries func (gs *GeospatialService) findGeographicalFeaturesWithinRadius(ctx context.Context, featureType domain.GeographicalFeatureType, lat, lng, radiusKm float64) ([]*domain.GeographicalFeature, error) { var features []*domain.GeographicalFeature query := ` SELECT * FROM geographical_features WHERE feature_type = ? AND ST_DWithin( geometry::geography, ST_SetSRID(ST_MakePoint(?, ?), 4326)::geography, ? * 1000 ) ORDER BY ST_Distance(geometry::geography, ST_SetSRID(ST_MakePoint(?::double precision, ?::double precision), 4326)::geography) ` result := gs.db.WithContext(ctx).Raw(query, featureType, lng, lat, radiusKm, lng, lat).Scan(&features) if result.Error != nil { return nil, result.Error } return features, nil } // CalculateSiteEnvironmentalScore calculates environmental score based on nearby green spaces func (gs *GeospatialService) CalculateSiteEnvironmentalScore(ctx context.Context, siteLat, siteLng float64) (float64, error) { // If no geographical features repo is available (e.g. lightweight tests), // return a neutral default score rather than causing a nil-pointer panic. if gs.geoFeatureRepo == nil { return 5.0, nil } // Get green spaces within 2km greenSpaces, err := gs.geoFeatureRepo.GetGreenSpacesWithinRadius(ctx, siteLat, siteLng, 2.0) if err != nil { return 0, err } // Calculate score based on proximity and size of green spaces var totalScore float64 for range greenSpaces { // Calculate distance to green space (we'd need to add distance calculation) // For now, use a simple scoring based on count totalScore += 1.0 } // Normalize score (max 10 points for environmental rating) if totalScore > 10 { totalScore = 10 } return totalScore, nil } // CalculateTransportationAccessibility calculates accessibility score based on road network func (gs *GeospatialService) CalculateTransportationAccessibility(ctx context.Context, siteLat, siteLng float64) (float64, error) { // Get roads within 1km roads, err := gs.geoFeatureRepo.GetRoadsWithinRadius(ctx, siteLat, siteLng, 1.0) if err != nil { return 0, err } // Calculate accessibility based on road density roadCount := len(roads) // Simple scoring: more roads = better accessibility var score float64 switch { case roadCount >= 10: score = 10.0 case roadCount >= 5: score = 7.5 case roadCount >= 2: score = 5.0 case roadCount >= 1: score = 2.5 default: score = 0.0 } return score, nil } // GetGeographicalFeatureStatistics returns comprehensive statistics about geographical features func (gs *GeospatialService) GetGeographicalFeatureStatistics(ctx context.Context) (map[string]interface{}, error) { stats := make(map[string]interface{}) // Get counts by feature type featureTypes := []domain.GeographicalFeatureType{ domain.GeographicalFeatureTypeRoad, domain.GeographicalFeatureTypeGreenSpace, domain.GeographicalFeatureTypePOI, domain.GeographicalFeatureTypeRailway, domain.GeographicalFeatureTypeWater, domain.GeographicalFeatureTypeLandUse, } for _, featureType := range featureTypes { features, err := gs.geoFeatureRepo.GetByType(ctx, featureType) if err == nil { stats[string(featureType)] = map[string]interface{}{ "count": len(features), } } } // Get road network statistics roadStats, err := gs.geoFeatureRepo.GetRoadNetworkStatistics(ctx) if err == nil { stats["road_network"] = roadStats } return stats, nil } // FindOptimalFacilityLocations finds optimal locations for new facilities based on criteria func (gs *GeospatialService) FindOptimalFacilityLocations(ctx context.Context, criteria FacilityLocationCriteria) ([]FacilityLocation, error) { // Algorithm considers: // - Proximity to existing industrial sites // - Access to road network // - Distance from residential areas // - Environmental constraints // - Available utilities var locations []FacilityLocation // Query for areas with good road access and proximity to existing sites query := ` WITH candidate_areas AS ( SELECT ST_Buffer(s.location_geometry, 1000) as area, s.id as nearby_site_id, s.latitude, s.longitude FROM sites s WHERE s.location_geometry IS NOT NULL LIMIT 10 ) SELECT ST_AsText(ST_Centroid(area)) as center_point, COUNT(*) as nearby_sites, ST_Y(ST_Centroid(area)) as lat, ST_X(ST_Centroid(area)) as lng FROM candidate_areas GROUP BY area HAVING COUNT(*) >= ? LIMIT ? ` rows, err := gs.db.WithContext(ctx).Raw(query, criteria.MinNearbySites, criteria.MaxResults).Rows() if err != nil { return nil, err } defer rows.Close() for rows.Next() { var loc FacilityLocation var centerPoint string err := rows.Scan(¢erPoint, &loc.NearbySites, &loc.Latitude, &loc.Longitude) if err != nil { continue } // Calculate scores for this location envScore, _ := gs.CalculateSiteEnvironmentalScore(ctx, loc.Latitude, loc.Longitude) transportScore, _ := gs.CalculateTransportationAccessibility(ctx, loc.Latitude, loc.Longitude) loc.EnvironmentalScore = envScore loc.TransportationScore = transportScore loc.OverallScore = (envScore + transportScore) / 2.0 locations = append(locations, loc) } return locations, nil } // FacilityLocationCriteria defines criteria for optimal facility location search type FacilityLocationCriteria struct { MinNearbySites int `json:"min_nearby_sites"` MaxDistanceKm float64 `json:"max_distance_km"` RequireRoadAccess bool `json:"require_road_access"` MinEnvironmentalScore float64 `json:"min_environmental_score"` MaxResults int `json:"max_results"` } // FacilityLocation represents a potential facility location with scores type FacilityLocation struct { Latitude float64 `json:"latitude"` Longitude float64 `json:"longitude"` NearbySites int `json:"nearby_sites"` EnvironmentalScore float64 `json:"environmental_score"` TransportationScore float64 `json:"transportation_score"` OverallScore float64 `json:"overall_score"` }