import math
from qgis.core import (QgsCoordinateReferenceSystem, QgsApplication, 
                       QgsCoordinateReferenceSystemRegistry, QgsRectangle)

class CRSAnalyzer:
    def __init__(self):
        self.crs_registry = QgsApplication.coordinateReferenceSystemRegistry()
        self.crs_db_records = None
        self._load_crs_database()
    
    def _load_crs_database(self):
        """Load all CRS records from QGIS internal database"""
        try:
            self.crs_db_records = self.crs_registry.crsDbRecords()[1]
        except:
            self.crs_db_records = []
    
    def detect_coordinate_type(self, extent):
        """
        Detect if coordinates are in decimal degrees or meters based on extent values
        
        Args:
            extent: QgsRectangle with layer extent
            
        Returns:
            dict with coordinate type information
        """
        if not extent or extent.isEmpty():
            return {'type': 'unknown', 'confidence': 0, 'reason': 'Invalid extent'}
        
        min_x, min_y = extent.xMinimum(), extent.yMinimum()
        max_x, max_y = extent.xMaximum(), extent.yMaximum()
        
        # Calculate center and ranges
        center_x = (min_x + max_x) / 2
        center_y = (min_y + max_y) / 2
        range_x = abs(max_x - min_x)
        range_y = abs(max_y - min_y)
        
        # Geographic coordinate detection (decimal degrees)
        if (abs(min_x) <= 180 and abs(max_x) <= 180 and 
            abs(min_y) <= 90 and abs(max_y) <= 90):
            
            # Additional checks for geographic coordinates
            if (range_x <= 360 and range_y <= 180):  # Within global bounds
                confidence = 95
                
                # Higher confidence if values look like typical lat/lon
                if (abs(center_x) <= 180 and abs(center_y) <= 90 and
                    0.001 <= range_x <= 360 and 0.001 <= range_y <= 180):
                    confidence = 98
                
                return {
                    'type': 'geographic',
                    'units': 'degrees',
                    'confidence': confidence,
                    'reason': f'Coordinates within geographic bounds (±180,±90)',
                    'center': (center_x, center_y),
                    'range': (range_x, range_y)
                }
        
        # Projected coordinate detection (meters/feet)
        # Check for large coordinate values typical of projected systems
        if (abs(min_x) > 180 or abs(max_x) > 180 or 
            abs(min_y) > 90 or abs(max_y) > 90):
            
            # UTM-like coordinates (6-7 digit eastings, 7-8 digit northings)
            if (10000 <= abs(center_x) <= 1000000 and 
                100000 <= abs(center_y) <= 10000000):
                return {
                    'type': 'projected',
                    'units': 'meters',
                    'confidence': 92,
                    'reason': 'UTM-like coordinate ranges',
                    'center': (center_x, center_y),
                    'range': (range_x, range_y),
                    'likely_system': 'UTM'
                }
            
            # Large projected coordinates (like Indian examples)
            elif (abs(center_x) > 1000000 or abs(center_y) > 1000000):
                return {
                    'type': 'projected',
                    'units': 'meters',
                    'confidence': 90,
                    'reason': 'Large coordinate values suggest projected system',
                    'center': (center_x, center_y),
                    'range': (range_x, range_y),
                    'likely_system': 'National/Regional Grid'
                }
            
            # Medium range projected coordinates
            elif (1000 <= abs(center_x) <= 10000000 and 
                  1000 <= abs(center_y) <= 10000000):
                return {
                    'type': 'projected',
                    'units': 'meters',
                    'confidence': 85,
                    'reason': 'Coordinate magnitude suggests projected system',
                    'center': (center_x, center_y),
                    'range': (range_x, range_y)
                }
            
            # Very large coordinates (could be feet)
            elif (abs(center_x) > 10000000 or abs(center_y) > 10000000):
                return {
                    'type': 'projected',
                    'units': 'feet',
                    'confidence': 80,
                    'reason': 'Very large coordinates suggest feet-based system',
                    'center': (center_x, center_y),
                    'range': (range_x, range_y)
                }
        
        # Small coordinates that might be in a local system
        if (abs(center_x) < 1000 and abs(center_y) < 1000 and
            range_x < 1000 and range_y < 1000):
            return {
                'type': 'local',
                'units': 'meters',
                'confidence': 70,
                'reason': 'Small coordinates suggest local coordinate system',
                'center': (center_x, center_y),
                'range': (range_x, range_y)
            }
        
        # Default case - uncertain
        return {
            'type': 'unknown',
            'units': 'unknown',
            'confidence': 30,
            'reason': 'Could not determine coordinate type from extent',
            'center': (center_x, center_y),
            'range': (range_x, range_y)
        }
    
    def get_suitable_crs(self, extent, is_vector=True, max_results=15):
        """
        Analyze layer extent and return suitable CRS options with intelligent coordinate type detection
        """
        if not extent:
            return self.get_common_crs()
        
        # Detect coordinate type from extent
        coord_info = self.detect_coordinate_type(extent)
        
        # Get all CRS from database
        all_crs = self._get_all_crs_from_database()
        suitable_crs = []
        
        # Calculate extent properties for geographic analysis
        width = extent.width()
        height = extent.height()
        center_lon = extent.center().x()
        center_lat = extent.center().y()
        area = width * height
        extent_type = self._classify_extent(area, coord_info['type'])
        
        for crs_record in all_crs:
            if self._is_crs_suitable_for_extent(crs_record, extent, extent_type, coord_info):
                crs_info = self._convert_db_record_to_info(crs_record)
                if crs_info:
                    suitable_crs.append(crs_info)
        
        # Score CRS with coordinate type intelligence
        scored_crs = []
        for crs_info in suitable_crs:
            score = self._calculate_intelligent_crs_score(
                crs_info, extent, is_vector, extent_type, coord_info)
            crs_info['score'] = score
            crs_info['coord_detection'] = coord_info
            scored_crs.append(crs_info)
        
        # Sort by score and return top results
        scored_crs.sort(key=lambda x: x['score'], reverse=True)
        return scored_crs[:max_results]
    
    def _classify_extent(self, area, coord_type):
        """Classify extent based on area and coordinate type"""
        if coord_type == 'geographic':
            # For geographic coordinates (degrees)
            if area > 14400:  # > 120x120 degrees
                return 'world'
            elif area > 900:  # > 30x30 degrees
                return 'continental'
            else:
                return 'regional'
        else:
            # For projected coordinates (meters/feet)
            if area > 1e14:  # Very large area in square meters
                return 'world'
            elif area > 1e12:  # Large area
                return 'continental'
            else:
                return 'regional'
    
    def _calculate_intelligent_crs_score(self, crs_info, extent, is_vector, extent_type, coord_info):
        """Calculate CRS score with intelligent coordinate type matching"""
        score = 30  # Base score (reduced from 50)
        
        auth_id = crs_info['authid']
        description = crs_info['description'].lower()
        crs_type = crs_info['type']
        coord_type = coord_info['type']
        coord_confidence = coord_info['confidence']
        
        # MAJOR BONUS: Coordinate type matching
        if coord_type == 'geographic' and crs_type == 'geographic':
            # Geographic coordinates with geographic CRS
            score += 80  # Major bonus
        elif coord_type in ['projected', 'local'] and crs_type == 'projected':
            # Projected coordinates with projected CRS
            score += 80  # Major bonus
        elif coord_type == 'geographic' and crs_type == 'projected':
            # Geographic coordinates but projected CRS suggested - penalty
            score -= 40
        elif coord_type in ['projected', 'local'] and crs_type == 'geographic':
            # Projected coordinates but geographic CRS suggested - major penalty
            score -= 60
        
        # Confidence bonus
        score += (coord_confidence / 10)  # Up to 10 bonus points
        
        # Unit-specific bonuses
        if coord_info.get('units') == 'meters':
            if 'meter' in description or 'metre' in description:
                score += 25
            elif 'feet' in description or 'foot' in description:
                score -= 15  # Penalty for unit mismatch
        
        # UTM-specific logic for projected coordinates
        if coord_type in ['projected', 'local'] and 'utm' in description:
            if coord_info.get('likely_system') == 'UTM':
                score += 40  # High bonus for UTM detection
            
            # Calculate optimal UTM zone if geographic center can be estimated
            center_x, center_y = coord_info['center']
            
            # Try to estimate geographic center for UTM zone calculation
            if coord_type == 'projected' and coord_info.get('likely_system') == 'UTM':
                # Rough UTM zone estimation from easting
                if 100000 <= abs(center_x) <= 1000000:
                    estimated_zone = max(1, min(60, int(abs(center_x) / 100000)))
                    if self._is_utm_zone_match(auth_id, estimated_zone):
                        score += 30
        
        # Geographic extent analysis (if coordinates allow)
        if coord_type == 'geographic':
            center_lon = extent.center().x()
            center_lat = extent.center().y()
            
            # Regional CRS bonuses for geographic coordinates
            if self._is_regional_crs_suitable(auth_id, description, center_lon, center_lat):
                score += 35
            
            # UTM zone calculation for geographic coordinates
            if 'utm' in description:
                optimal_zone = int((center_lon + 180) / 6) + 1
                if self._is_utm_zone_match(auth_id, optimal_zone):
                    score += 50  # Perfect UTM zone
        
        # Common CRS handling
        if auth_id == 'EPSG:4326':  # WGS84
            if coord_type == 'geographic':
                score += 40  # High bonus for geographic data
            else:
                score -= 20  # Penalty for projected data
        elif auth_id == 'EPSG:3857':  # Web Mercator
            if coord_type == 'projected':
                score += 30
            else:
                score -= 10
        
        # Extent type bonuses
        if extent_type == 'world':
            world_keywords = ['world', 'global', 'mollweide', 'robinson']
            if any(keyword in description for keyword in world_keywords):
                score += 30
        elif extent_type == 'regional':
            if 'utm' in description or 'state plane' in description:
                score += 25
        
        # Deprecated CRS penalty
        deprecated_keywords = ['nad27', 'ed50', 'deprecated']
        if any(keyword in description for keyword in deprecated_keywords):
            score -= 25
        
        return max(0, score)
    
    def _is_utm_zone_match(self, auth_id, zone):
        """Check if EPSG code matches the expected UTM zone"""
        try:
            if auth_id.startswith('EPSG:326') and len(auth_id) == 10:  # Northern hemisphere
                crs_zone = int(auth_id[-2:])
                return abs(crs_zone - zone) <= 1  # Allow neighboring zones
            elif auth_id.startswith('EPSG:327') and len(auth_id) == 10:  # Southern hemisphere
                crs_zone = int(auth_id[-2:])
                return abs(crs_zone - zone) <= 1
        except:
            pass
        return False
    
    def _get_all_crs_from_database(self):
        """Get all CRS records from QGIS database"""
        all_crs = []
        
        # Get common EPSG codes with expanded coverage
        epsg_codes = self._get_comprehensive_epsg_codes()
        
        for code in epsg_codes:
            try:
                crs = QgsCoordinateReferenceSystem(f"EPSG:{code}")
                if crs.isValid():
                    crs_record = {
                        'auth_id': f"EPSG:{code}",
                        'description': crs.description(),
                        'is_geographic': crs.isGeographic(),
                        'crs': crs
                    }
                    all_crs.append(crs_record)
            except:
                continue
        
        return all_crs
    
    def _get_comprehensive_epsg_codes(self):
        """Comprehensive list of EPSG codes"""
        # Geographic CRS (major datums)
        geographic_codes = [4326, 4269, 4258, 4167, 4283, 4314, 4019, 4218, 4230, 4267]
        
        # Major projected CRS
        projected_codes = [3857, 3035, 3031, 2154, 5070, 3413, 3574, 3577, 27700, 2163]
        
        # UTM zones - complete coverage
        utm_north = list(range(32601, 32661))  # All northern UTM zones
        utm_south = list(range(32701, 32761))  # All southern UTM zones
        
        # Regional systems by continent/country
        
        # North America
        na_codes = list(range(2225, 2290))  # State Plane NAD83
        na_codes.extend(range(3311, 3370))  # Additional state plane
        na_codes.extend([5070, 5072, 2163])  # Albers and other NA systems
        
        # Europe
        eu_codes = list(range(25828, 25838))  # ETRS89 UTM zones
        eu_codes.extend(range(3038, 3051))    # ETRS89 TM zones
        eu_codes.extend([3035, 3857, 4258])   # ETRS89 LAEA, Web Mercator
        
        # United Kingdom and Ireland
        uk_codes = [27700, 29901, 29902, 7405]  # BNG, Irish Grid
        
        # Australia and New Zealand
        au_codes = list(range(28348, 28358))  # GDA94 MGA zones
        au_codes.extend([3577, 7844, 7845, 2193])  # Albers, NZTM
        
        # Asia
        asia_codes = list(range(3414, 3440))  # Various Asian systems
        asia_codes.extend([3857, 4326])        # Common global systems
        
        # Africa
        africa_codes = list(range(32628, 32638))  # UTM zones for Africa
        africa_codes.extend([4326, 3857])         # Global systems
        
        # South America
        sa_codes = list(range(32717, 32725))  # UTM zones for South America
        sa_codes.extend([4326, 3857])         # Global systems
        
        # India specific (based on user's example)
        india_codes = [32643, 32644, 32645, 32646, 7755, 7756]  # UTM and local systems
        
        # Combine all codes
        all_codes = (geographic_codes + projected_codes + utm_north + utm_south + 
                    na_codes + eu_codes + uk_codes + au_codes + asia_codes + 
                    africa_codes + sa_codes + india_codes)
        
        return sorted(set(all_codes))
    
    def _convert_db_record_to_info(self, crs_record):
        """Convert database record to CRS info dictionary"""
        try:
            if isinstance(crs_record, dict):
                auth_id = crs_record.get('auth_id', '')
                description = crs_record.get('description', '')
                crs = crs_record.get('crs')
                
                if not crs:
                    crs = QgsCoordinateReferenceSystem(auth_id)
            else:
                auth_id = getattr(crs_record, 'authId', lambda: '')()
                description = getattr(crs_record, 'description', lambda: '')()
                crs = QgsCoordinateReferenceSystem(auth_id)
            
            if not crs.isValid():
                return None
            
            crs_type = 'geographic' if crs.isGeographic() else 'projected'
            
            return {
                'authid': auth_id,
                'description': description,
                'type': crs_type,
                'crs_object': crs
            }
            
        except Exception as e:
            return None
    
    def _is_crs_suitable_for_extent(self, crs_record, extent, extent_type, coord_info):
        """Enhanced CRS suitability check considering coordinate type"""
        try:
            crs_info = self._convert_db_record_to_info(crs_record)
            if not crs_info:
                return False
            
            auth_id = crs_info['authid']
            description = crs_info['description'].lower()
            crs_type = crs_info['type']
            coord_type = coord_info['type']
            
            # Always include WGS84 and Web Mercator
            if auth_id in ['EPSG:4326', 'EPSG:3857']:
                return True
            
            # Strong preference for matching coordinate types
            if coord_type == 'geographic' and crs_type != 'geographic':
                # Only include projected CRS if they're commonly used with geographic data
                if not any(keyword in description for keyword in ['utm', 'web mercator', 'pseudo']):
                    return False
            
            if coord_type in ['projected', 'local'] and crs_type != 'projected':
                # Rarely include geographic CRS for projected data
                return auth_id == 'EPSG:4326'  # Only WGS84 as fallback
            
            # Rest of the filtering logic...
            center_lon = extent.center().x()
            center_lat = extent.center().y()
            
            # For geographic coordinates, use geographic center directly
            if coord_type == 'geographic':
                if extent_type == 'world':
                    world_keywords = ['world', 'global', 'mollweide', 'robinson', 'eckert']
                    return any(keyword in description for keyword in world_keywords)
                
                elif extent_type == 'regional':
                    if 'utm' in description:
                        return self._is_utm_suitable(auth_id, center_lon, center_lat, buffer_zones=2)
                    
                    if self._is_regional_crs_suitable(auth_id, description, center_lon, center_lat):
                        return True
            
            # For projected coordinates, be more liberal with CRS selection
            else:
                if 'utm' in description:
                    return True  # Include all UTM zones for projected data
                
                # Include regional systems
                regional_keywords = ['albers', 'lambert', 'transverse mercator', 'state plane']
                if any(keyword in description for keyword in regional_keywords):
                    return True
            
            return False
            
        except:
            return False
    
    def _is_utm_suitable(self, auth_id, center_lon, center_lat, buffer_zones=1):
        """Check if UTM zone is suitable for the location"""
        try:
            optimal_zone = int((center_lon + 180) / 6) + 1
            
            if auth_id.startswith('EPSG:326'):  # Northern hemisphere
                zone = int(auth_id[-2:])
                hemisphere = 'N'
            elif auth_id.startswith('EPSG:327'):  # Southern hemisphere
                zone = int(auth_id[-2:])
                hemisphere = 'S'
            else:
                return False
            
            if (center_lat >= 0 and hemisphere != 'N') or (center_lat < 0 and hemisphere != 'S'):
                return False
            
            zone_diff = abs(zone - optimal_zone)
            return zone_diff <= buffer_zones
            
        except:
            return False
    
    def _is_regional_crs_suitable(self, auth_id, description, center_lon, center_lat):
        """Check if regional CRS is suitable based on geographic location"""
        # [Same implementation as before]
        if -130 <= center_lon <= -60 and 20 <= center_lat <= 70:
            na_keywords = ['nad83', 'nad27', 'conus', 'albers', 'state plane', 'spcs']
            if any(keyword in description for keyword in na_keywords):
                return True
        
        if -15 <= center_lon <= 45 and 35 <= center_lat <= 75:
            eu_keywords = ['etrs89', 'ed50', 'europe', 'laea', 'tmerc']
            if any(keyword in description for keyword in eu_keywords):
                return True
        
        if 110 <= center_lon <= 180 and -50 <= center_lat <= -5:
            au_keywords = ['gda94', 'gda2020', 'australian', 'albers', 'mga']
            if any(keyword in description for keyword in au_keywords):
                return True
        
        return False
    
    def get_common_crs(self):
        """Get common CRS options when extent analysis fails"""
        common_epsg_codes = [
            (4326, 'WGS 84 (Geographic)', 'geographic', 100),
            (3857, 'WGS 84 / Pseudo-Mercator (Web Mercator)', 'projected', 90),
            (4269, 'NAD83 (Geographic)', 'geographic', 80),
            (3035, 'ETRS89 / LAEA Europe', 'projected', 75),
            (5070, 'NAD83 / Conus Albers (USA)', 'projected', 75),
            (3577, 'GDA94 / Australian Albers', 'projected', 75),
            (32633, 'WGS 84 / UTM zone 33N', 'projected', 70),
            (32643, 'WGS 84 / UTM zone 43N (India)', 'projected', 70),
        ]
        
        common_crs = []
        for epsg_code, desc, crs_type, score in common_epsg_codes:
            crs = QgsCoordinateReferenceSystem(f'EPSG:{epsg_code}')
            if crs.isValid():
                common_crs.append({
                    'authid': f'EPSG:{epsg_code}',
                    'description': desc,
                    'type': crs_type,
                    'score': score,
                    'crs_object': crs
                })
        
        return common_crs

    def detect_likely_crs_from_extent(self, extent):
        """Detect likely CRS based on coordinate ranges"""
        min_x, min_y = extent.xMinimum(), extent.yMinimum()
        max_x, max_y = extent.xMaximum(), extent.yMaximum()
        
        likely_crs = []
        
        # Geographic coordinates (lat/lon)
        if (-180 <= min_x <= 180 and -180 <= max_x <= 180 and 
            -90 <= min_y <= 90 and -90 <= max_y <= 90):
            crs = QgsCoordinateReferenceSystem('EPSG:4326')
            likely_crs.append({
                'authid': 'EPSG:4326', 
                'description': 'WGS 84 (Geographic)', 
                'type': 'geographic',
                'score': 95,
                'reason': 'Coordinates within geographic range',
                'crs_object': crs
            })
        
        # UTM-like coordinates
        center_x = (min_x + max_x) / 2
        center_y = (min_y + max_y) / 2
        
        if (100000 <= abs(center_x) <= 1000000 and 
            1000000 <= abs(center_y) <= 10000000):
            
            # Estimate UTM zone from easting
            estimated_zone = max(1, min(60, int(abs(center_x) / 1000000) + 30))
            
            if center_y > 0:
                utm_epsg = 32600 + estimated_zone
                utm_desc = f'WGS 84 / UTM zone {estimated_zone}N'
            else:
                utm_epsg = 32700 + estimated_zone
                utm_desc = f'WGS 84 / UTM zone {estimated_zone}S'
            
            crs = QgsCoordinateReferenceSystem(f'EPSG:{utm_epsg}')
            if crs.isValid():
                likely_crs.append({
                    'authid': f'EPSG:{utm_epsg}',
                    'description': utm_desc,
                    'type': 'projected',
                    'score': 85,
                    'reason': 'Coordinates suggest UTM projection',
                    'crs_object': crs
                })
        
        return likely_crs
