nakas/NWPS_SWAN
0
1import os2import sys3import tempfile4import logging5import subprocess6import shutil7from datetime import datetime, timedelta8import numpy as np9import xarray as xr10from ecmwf.opendata import Client11import requests12 13# Setup logging first14logging.basicConfig(level=logging.INFO)15logger = logging.getLogger(__name__)16 17# Add current directory to path for Arctic extractor import18sys.path.append(os.path.dirname(os.path.abspath(__file__)))19 20# Import ONLY the working Arctic GRIB handler for Docker production21try:22 from arctic_grib_handler import ArcticGRIBHandler23 ARCTIC_HANDLER_AVAILABLE = True24 logger.info("✅ Arctic GRIB Handler loaded (Docker production version)")25except ImportError as e:26 logger.error(f"❌ Arctic GRIB handler not available: {e}")27 ARCTIC_HANDLER_AVAILABLE = False28 29class GRIBWavePuller:30 def __init__(self):31 self.client = Client("ecmwf")32 self.output_dir = os.getenv('OUTPUT_DIR', '/tmp/wave_data')33 os.makedirs(self.output_dir, exist_ok=True)34 35 # Set ECCODES environment variables to handle polar stereographic issues36 self._setup_eccodes_environment()37 38 def _setup_eccodes_environment(self):39 """Setup ECCODES environment variables to handle projection issues"""40 try:41 # Set environment variables that might help with polar stereographic processing42 os.environ['ECCODES_GRIB_STRICT_PARSING'] = '0' # Relaxed parsing43 os.environ['ECCODES_GRIB_IGNORE_GRID_DEFINITION'] = '1' # Ignore grid definition errors44 logger.info("Set ECCODES environment variables for relaxed parsing")45 except Exception as e:46 logger.warning(f"Could not set ECCODES environment variables: {e}")47 48 def _check_cdo_available(self):49 """Check if CDO (Climate Data Operators) is available"""50 try:51 result = subprocess.run(['cdo', '--version'], 52 capture_output=True, text=True, timeout=10)53 return result.returncode == 054 except (subprocess.TimeoutExpired, FileNotFoundError, subprocess.SubprocessError):55 return False56 57 def _reproject_arctic_with_cdo(self, grib_file_path):58 """Reproject Arctic GRIB file using CDO as alternative to wgrib2"""59 try:60 if not self._check_cdo_available():61 logger.warning("CDO not available for Arctic reprojection")62 return None63 64 logger.info("Attempting to reproject Arctic GRIB file using CDO")65 66 # Create temporary file for reprojected data67 temp_reprojected = tempfile.NamedTemporaryFile(delete=False, suffix='_cdo_reprojected.grib2')68 temp_reprojected.close()69 70 # Use CDO to reproject to regular lat-lon grid71 # remapbil = bilinear interpolation to regular lat-lon grid72 cmd = [73 'cdo', 'remapbil,r720x360', # 0.5° resolution global grid74 grib_file_path,75 temp_reprojected.name76 ]77 78 logger.info(f"Running CDO command: {' '.join(cmd)}")79 result = subprocess.run(cmd, capture_output=True, text=True, timeout=120)80 81 if result.returncode == 0:82 logger.info("Successfully reprojected Arctic GRIB file with CDO")83 return temp_reprojected.name84 else:85 logger.error(f"CDO failed: {result.stderr}")86 if os.path.exists(temp_reprojected.name):87 os.unlink(temp_reprojected.name)88 return None89 90 except Exception as e:91 logger.error(f"Error reprojecting with CDO: {e}")92 return None93 94 def fetch_ecmwf_wave_grib(self, forecast_time=0):95 """Fetch global wave data from ECMWF open data"""96 try:97 logger.info("Fetching ECMWF global wave GRIB data...")98 99 # Create temporary file for GRIB data100 temp_file = tempfile.NamedTemporaryFile(delete=False, suffix='.grib2')101 102 try:103 # Fetch wave height data from ECMWF104 self.client.retrieve(105 type="fc", # forecast106 param=["swh"], # significant wave height107 time=0, # 00 UTC108 step=forecast_time, # forecast hours ahead109 target=temp_file.name110 )111 112 logger.info(f"GRIB file downloaded: {temp_file.name}")113 return temp_file.name114 115 except Exception as e:116 logger.error(f"Failed to fetch ECMWF data: {e}")117 # Clean up temp file on error118 if os.path.exists(temp_file.name):119 os.unlink(temp_file.name)120 return None121 122 except Exception as e:123 logger.error(f"Error in fetch_ecmwf_wave_grib: {e}")124 return None125 126 def fetch_noaa_wave_grib(self, forecast_hour=0):127 """Fetch global wave data from NOAA WW3 model"""128 try:129 logger.info(f"Fetching NOAA WW3 global wave GRIB data for forecast hour {forecast_hour}...")130 131 # NOAA GFS/WW3 wave data URL pattern132 base_url = "https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod"133 134 # Try current date and previous days (in case of delayed updates)135 now = datetime.utcnow()136 dates_to_try = [137 now.strftime("%Y%m%d"),138 (now - timedelta(days=1)).strftime("%Y%m%d"),139 (now - timedelta(days=2)).strftime("%Y%m%d")140 ]141 142 # Try different model runs (00, 06, 12, 18 UTC) to find available data143 model_runs = ["00", "06", "12", "18"]144 current_hour = now.hour145 146 # Start with the most recent available run147 if current_hour >= 18:148 preferred_runs = ["18", "12", "06", "00"]149 elif current_hour >= 12:150 preferred_runs = ["12", "06", "00", "18"]151 elif current_hour >= 6:152 preferred_runs = ["06", "00", "18", "12"]153 else:154 preferred_runs = ["00", "18", "12", "06"]155 156 # Try different dates and model runs157 for date_str in dates_to_try:158 logger.info(f"Trying date: {date_str}")159 for hour in preferred_runs:160 try:161 # Format forecast hour with leading zeros (f000, f001, f002, etc.)162 forecast_str = f"f{forecast_hour:03d}"163 164 # Download multiple regional files for global coverage165 successful_downloads = []166 167 # Try different regional GRIB files available on NOAA168 regional_files = [169 (f"gfswave.t{hour}z.atlocn.0p16.{forecast_str}.grib2", "Atlantic"),170 (f"gfswave.t{hour}z.epacif.0p16.{forecast_str}.grib2", "East_Pacific"), 171 (f"gfswave.t{hour}z.arctic.9km.{forecast_str}.grib2", "Arctic"),172 # Add more regional files if available173 (f"gfswave.t{hour}z.global.0p16.{forecast_str}.grib2", "Global"), # Try global if it exists174 ]175 176 # Try to download each regional file177 for filename, region_name in regional_files:178 try:179 url = f"{base_url}/gfs.{date_str}/{hour}/wave/gridded/{filename}"180 logger.info(f"Attempting to download {region_name} region: {filename}")181 182 temp_file = tempfile.NamedTemporaryFile(delete=False, suffix='.grib2')183 184 response = requests.get(url, timeout=300)185 if response.status_code == 200:186 temp_file.write(response.content)187 temp_file.close()188 successful_downloads.append((temp_file.name, region_name, hour, forecast_hour))189 logger.info(f"{region_name} GRIB file downloaded: {temp_file.name}")190 else:191 logger.debug(f"HTTP {response.status_code} for {region_name}")192 os.unlink(temp_file.name)193 continue194 except Exception as file_error:195 logger.debug(f"Error downloading {region_name}: {file_error}")196 continue197 198 # If we got at least one regional file, return the list199 if successful_downloads:200 logger.info(f"Successfully downloaded {len(successful_downloads)} regional files")201 return successful_downloads202 203 except Exception as run_error:204 logger.warning(f"Error trying {hour}Z run on {date_str}: {run_error}")205 continue206 207 logger.error("Failed to download NOAA data from any model run")208 return None, None, None209 210 except Exception as e:211 logger.error(f"Error in fetch_noaa_wave_grib: {e}")212 return None, None, None213 214 def _check_wgrib2_available(self):215 """Check if wgrib2 command-line tool is available"""216 try:217 result = subprocess.run(['wgrib2', '-version'], 218 capture_output=True, text=True, timeout=10)219 return result.returncode == 0220 except (subprocess.TimeoutExpired, FileNotFoundError, subprocess.SubprocessError):221 return False222 223 def _reproject_arctic_with_wgrib2(self, grib_file_path):224 """Reproject Arctic GRIB file to lat-lon grid using wgrib2"""225 try:226 if not self._check_wgrib2_available():227 logger.warning("wgrib2 not available for Arctic reprojection")228 return None229 230 logger.info("Attempting to reproject Arctic GRIB file using wgrib2")231 232 # Create temporary file for reprojected data233 temp_reprojected = tempfile.NamedTemporaryFile(delete=False, suffix='_reprojected.grib2')234 temp_reprojected.close()235 236 # Use wgrib2 to reproject to lat-lon grid237 # This covers Arctic regions with reasonable resolution238 cmd = [239 'wgrib2', grib_file_path,240 '-new_grid', 'latlon', '0:720:0.5', '50:71:0.5', # 0.5° res, 50-85°N, 0-360°E241 temp_reprojected.name242 ]243 244 logger.info(f"Running wgrib2 command: {' '.join(cmd)}")245 result = subprocess.run(cmd, capture_output=True, text=True, timeout=120)246 247 if result.returncode == 0:248 logger.info("Successfully reprojected Arctic GRIB file")249 return temp_reprojected.name250 else:251 logger.error(f"wgrib2 failed: {result.stderr}")252 if os.path.exists(temp_reprojected.name):253 os.unlink(temp_reprojected.name)254 return None255 256 except Exception as e:257 logger.error(f"Error reprojecting with wgrib2: {e}")258 return None259 260 261 def _process_arctic_without_coordinates(self, grib_file_path):262 """Process Arctic GRIB file bypassing coordinate processing entirely"""263 try:264 logger.info("Attempting to process Arctic file without coordinate processing")265 266 # Try to open with minimal processing - just get the data values267 try:268 # Use cfgrib with very restrictive read_keys to bypass coordinate issues269 ds = xr.open_dataset(270 grib_file_path, 271 engine='cfgrib',272 decode_timedelta=True,273 backend_kwargs={274 'read_keys': ['paramId', 'shortName', 'name', 'units'], # Only read essential keys275 'errors': 'ignore',276 'indexpath': ''277 }278 )279 280 logger.info(f"Successfully opened Arctic file with restricted processing")281 logger.info(f"Available variables: {list(ds.variables.keys())}")282 283 # Look for wave height data284 wave_var = None285 for var_name in ['swh', 'HTSGW', 'htsgw']:286 if var_name in ds.variables:287 wave_var = var_name288 break289 290 if wave_var is None:291 # Try broader search292 for var_name in ds.variables:293 if any(keyword in var_name.lower() for keyword in ['wave', 'height', 'swh']):294 wave_var = var_name295 break296 297 if wave_var:298 wave_data = ds[wave_var].values299 logger.info(f"Found wave data: {wave_var}, shape: {wave_data.shape}")300 301 # Create fake coordinate grid for Arctic region (approximate)302 # This is a fallback when we can't get real coordinates303 if len(wave_data.shape) == 2:304 rows, cols = wave_data.shape305 # Create approximate Arctic coordinate grid306 fake_lats = np.linspace(85, 60, rows) # 85°N to 60°N307 fake_lons = np.linspace(-180, 180, cols) # Full longitude range308 309 lon_grid, lat_grid = np.meshgrid(fake_lons, fake_lats)310 311 # Flatten and filter valid data312 flat_lats = lat_grid.flatten()313 flat_lons = lon_grid.flatten()314 flat_waves = wave_data.flatten()315 316 # Remove invalid data317 valid_mask = (~np.isnan(flat_waves)) & (flat_waves >= 0) & (flat_waves < 50)318 319 if np.any(valid_mask):320 filtered_lats = flat_lats[valid_mask]321 filtered_lons = flat_lons[valid_mask]322 filtered_waves = flat_waves[valid_mask]323 324 logger.info(f"Arctic fallback: {len(filtered_lats)} approximate points")325 return filtered_lats, filtered_lons, filtered_waves, None, None326 327 ds.close()328 return None, None329 330 except Exception as restricted_error:331 logger.warning(f"Restricted processing failed: {restricted_error}")332 return None, None333 334 except Exception as e:335 logger.error(f"Error in coordinate-bypass processing: {e}")336 return None, None337 338 def process_grib_file(self, grib_file_path, region_name=None):339 """Process GRIB file and extract wave data including direction and period"""340 try:341 logger.info(f"Processing GRIB file: {grib_file_path}")342 343 # Check if this is an Arctic file that needs special handling344 is_arctic = (region_name and 'arctic' in region_name.lower()) or 'arctic' in grib_file_path.lower()345 346 # Try to open GRIB file and extract all available wave parameters347 try:348 datasets = []349 350 if is_arctic:351 # Multi-layered approach for Arctic polar stereographic projection352 logger.info("Processing Arctic GRIB file with enhanced multi-layered fallback approach")353 354 # Approach 1: Try wgrib2 reprojection first (most reliable)355 reprojected_file = self._reproject_arctic_with_wgrib2(grib_file_path)356 if reprojected_file:357 try:358 logger.info("Processing wgrib2-reprojected Arctic file")359 ds_height = xr.open_dataset(reprojected_file, engine='cfgrib',360 decode_timedelta=True)361 datasets.append(ds_height)362 # Clean up reprojected file after processing363 os.unlink(reprojected_file)364 except Exception as reprojected_error:365 logger.warning(f"Failed to process reprojected Arctic file: {reprojected_error}")366 if os.path.exists(reprojected_file):367 os.unlink(reprojected_file)368 datasets = []369 370 # Approach 2: If wgrib2 failed, try CDO reprojection371 if not datasets:372 cdo_reprojected_file = self._reproject_arctic_with_cdo(grib_file_path)373 if cdo_reprojected_file:374 try:375 logger.info("Processing CDO-reprojected Arctic file")376 ds_height = xr.open_dataset(cdo_reprojected_file, engine='cfgrib',377 decode_timedelta=True)378 datasets.append(ds_height)379 # Clean up reprojected file after processing380 os.unlink(cdo_reprojected_file)381 except Exception as cdo_error:382 logger.warning(f"Failed to process CDO-reprojected Arctic file: {cdo_error}")383 if os.path.exists(cdo_reprojected_file):384 os.unlink(cdo_reprojected_file)385 datasets = []386 387 # Approach 3: Try coordinate-bypass method388 if not datasets:389 try:390 logger.info("Trying coordinate-bypass method for Arctic processing")391 result = self._process_arctic_without_coordinates(grib_file_path)392 if result and result[0] is not None:393 return result394 except Exception as bypass_error:395 logger.warning(f"Coordinate-bypass method failed: {bypass_error}")396 397 # Approach 4: If all reprojection methods failed, try cfgrib with relaxed settings398 if not datasets:399 try:400 logger.info("Trying cfgrib with relaxed error handling for Arctic")401 ds_height = xr.open_dataset(grib_file_path, engine='cfgrib',402 backend_kwargs={'errors': 'ignore'},403 decode_timedelta=True)404 datasets.append(ds_height)405 logger.info("Successfully opened Arctic file with relaxed cfgrib settings")406 except Exception as cfgrib_error:407 logger.warning(f"Failed to process Arctic file with cfgrib: {cfgrib_error}")408 409 # Use ONLY the proven working Arctic GRIB handler410 if not datasets:411 try:412 logger.info("🧊 Using proven Arctic GRIB handler (Docker production)")413 414 if not ARCTIC_HANDLER_AVAILABLE:415 logger.error("❌ Arctic GRIB handler not available")416 return None, None417 418 arctic_handler = ArcticGRIBHandler()419 result = arctic_handler.get_compatible_format(grib_file_path, sample_points=200)420 421 if result['status'] == 'success' and result['sampled_points'] > 0:422 logger.info(f"✅ Arctic extraction successful: {result['sampled_points']} points")423 data = result['data']424 # Convert to expected format (latitude, longitude, value)425 return (data['latitude'], data['longitude'], data['value'], None, None)426 else:427 logger.error(f"❌ Arctic extraction failed: {result['message']}")428 return None, None429 430 except Exception as extraction_error:431 logger.error(f"❌ Arctic GRIB handler failed: {extraction_error}")432 return None, None433 else:434 # Normal processing for non-Arctic files with timedelta fix435 ds_height = xr.open_dataset(grib_file_path, engine='cfgrib',436 decode_timedelta=True)437 datasets.append(ds_height)438 439 # Try to get wave direction and period by opening with different filters440 try:441 ds_ocean = xr.open_dataset(grib_file_path, engine='cfgrib', 442 filter_by_keys={'discipline': 10},443 decode_timedelta=True)444 if ds_ocean.variables.keys() != ds_height.variables.keys():445 datasets.append(ds_ocean)446 except:447 logger.info("Could not open oceanographic discipline data separately")448 449 # Combine all available variables450 all_vars = {}451 for ds in datasets:452 all_vars.update(ds.variables)453 454 logger.info(f"Available variables: {list(all_vars.keys())}")455 456 except Exception as e:457 error_msg = str(e)458 logger.error(f"Error opening GRIB file: {error_msg}")459 460 # Check if this is an Arctic polar stereographic error461 if is_arctic and any(keyword in error_msg.lower() for keyword in [462 'polar stereographic', 'spherical earth', 'geoiterator', 463 'geographic attributes', 'unable to create iterator'464 ]):465 logger.info("Detected Arctic polar stereographic error - using eccodes-based Arctic extraction")466 467 try:468 logger.info("🧊 Arctic error detected - using proven handler")469 470 if not ARCTIC_HANDLER_AVAILABLE:471 logger.error("❌ Arctic GRIB handler not available")472 return None, None473 474 arctic_handler = ArcticGRIBHandler()475 result = arctic_handler.get_compatible_format(grib_file_path, sample_points=200)476 477 if result['status'] == 'success' and result['sampled_points'] > 0:478 logger.info(f"✅ Arctic error handling successful: {result['sampled_points']} points")479 data = result['data']480 # Convert to expected format (latitude, longitude, value) 481 return (data['latitude'], data['longitude'], data['value'], None, None)482 else:483 logger.error(f"❌ Arctic error handling failed: {result['message']}")484 return None, None485 486 except Exception as arctic_error:487 logger.error(f"❌ Arctic handler error processing failed: {arctic_error}")488 return None, None489 else:490 # Non-Arctic error or different error type491 return None, None492 493 # Extract wave height data494 wave_height_var = None495 wave_heights = None496 for var_name in ['swh', 'HTSGW', 'htsgw']:497 if var_name in all_vars:498 wave_height_var = var_name499 wave_heights = all_vars[var_name].values500 logger.info(f"Using wave height variable: {wave_height_var}")501 break502 503 if wave_heights is None:504 # Try broader search505 for var_name in all_vars:506 if any(keyword in var_name.lower() for keyword in ['wave', 'height', 'swh']):507 wave_height_var = var_name508 wave_heights = all_vars[var_name].values509 logger.info(f"Found wave height variable: {wave_height_var}")510 break511 512 if wave_heights is None:513 logger.error("No wave height variables found in GRIB file")514 for ds in datasets:515 ds.close()516 return None, None517 518 # Extract wave direction data519 wave_directions = None520 wave_dir_var = None521 for var_name in ['dirpw', 'DIRPW', 'dp', 'wvdir', 'WVDIR', 'dir']:522 if var_name in all_vars:523 wave_dir_var = var_name524 wave_directions = all_vars[var_name].values525 logger.info(f"Found wave direction variable: {wave_dir_var}")526 break527 528 # Extract wave period data529 wave_periods = None530 wave_period_var = None531 for var_name in ['perpw', 'PERPW', 'tp', 'wvper', 'WVPER', 'per']:532 if var_name in all_vars:533 wave_period_var = var_name534 wave_periods = all_vars[var_name].values535 logger.info(f"Found wave period variable: {wave_period_var}")536 break537 538 # Get coordinates from the first dataset539 ds_main = datasets[0]540 lats = ds_main.latitude.values if 'latitude' in ds_main else ds_main.lat.values541 lons = ds_main.longitude.values if 'longitude' in ds_main else ds_main.lon.values542 543 # Log what we found544 if wave_directions is not None:545 logger.info(f"Wave directions shape: {wave_directions.shape}, range: {np.nanmin(wave_directions):.1f}-{np.nanmax(wave_directions):.1f} degrees")546 if wave_periods is not None:547 logger.info(f"Wave periods shape: {wave_periods.shape}, range: {np.nanmin(wave_periods):.1f}-{np.nanmax(wave_periods):.1f} seconds")548 549 # Create structured data550 processed_data = {551 'timestamp': datetime.utcnow().isoformat(),552 'data_source': 'ECMWF_GRIB' if 'ecmwf' in grib_file_path.lower() else 'NOAA_GRIB',553 'parameters_found': {554 'wave_height': wave_height_var,555 'wave_direction': wave_dir_var,556 'wave_period': wave_period_var,557 'has_velocity_components': wave_dir_var is not None558 },559 'grid_info': {560 'lat_min': float(np.min(lats)),561 'lat_max': float(np.max(lats)),562 'lon_min': float(np.min(lons)),563 'lon_max': float(np.max(lons)),564 'lat_resolution': float(lats[1] - lats[0]) if len(lats) > 1 else None,565 'lon_resolution': float(lons[1] - lons[0]) if len(lons) > 1 else None,566 'grid_shape': wave_heights.shape567 },568 'wave_statistics': {569 'max_wave_height': float(np.nanmax(wave_heights)),570 'min_wave_height': float(np.nanmin(wave_heights)),571 'mean_wave_height': float(np.nanmean(wave_heights)),572 'std_wave_height': float(np.nanstd(wave_heights))573 },574 'sample_points': self._extract_sample_points_with_vectors(lats, lons, wave_heights, wave_directions, wave_periods)575 }576 577 # Add direction/period statistics if available578 if wave_directions is not None:579 processed_data['direction_statistics'] = {580 'mean_direction': float(np.nanmean(wave_directions)),581 'direction_std': float(np.nanstd(wave_directions))582 }583 584 if wave_periods is not None:585 processed_data['period_statistics'] = {586 'max_period': float(np.nanmax(wave_periods)),587 'min_period': float(np.nanmin(wave_periods)),588 'mean_period': float(np.nanmean(wave_periods))589 }590 591 # Close all datasets592 for ds in datasets:593 ds.close()594 595 return processed_data, grib_file_path596 597 except Exception as e:598 logger.error(f"Error processing GRIB file: {e}")599 return None, None600 601 def process_multiple_regional_files(self, regional_files):602 """Process multiple regional GRIB files and combine data for global coverage"""603 try:604 logger.info(f"Processing {len(regional_files)} regional GRIB files for global coverage...")605 606 combined_sample_points = []607 all_wave_heights = []608 all_wave_directions = []609 all_wave_periods = []610 611 global_lat_min = float('inf')612 global_lat_max = float('-inf')613 global_lon_min = float('inf') 614 global_lon_max = float('-inf')615 616 parameters_found = {617 'wave_height': None,618 'wave_direction': None,619 'wave_period': None,620 'has_velocity_components': False621 }622 623 regions_processed = []624 625 for grib_file_path, region_name, model_run, forecast_hour in regional_files:626 try:627 logger.info(f"Processing {region_name} region: {grib_file_path}")628 629 # Process this regional file630 result = self.process_grib_file(grib_file_path, region_name=region_name)631 632 # Handle different return formats (normal processing vs Arctic pygrib)633 if result is None or (isinstance(result, tuple) and len(result) == 2 and result[0] is None):634 logger.warning(f"Failed to process {region_name} region")635 continue636 637 # Check if this is Arctic data with raw coordinates (from pygrib)638 if (isinstance(result, tuple) and len(result) == 5 and 639 not isinstance(result[0], dict)):640 # Arctic data: (lats, lons, heights, directions, periods)641 lats, lons, heights, directions, periods = result642 logger.info(f"Processing Arctic raw coordinate data: {len(lats)} points")643 644 # Convert to sample points format645 regional_points = []646 for i in range(len(lats)):647 point = {648 'lat': float(lats[i]),649 'lon': float(lons[i]),650 'wave_height': float(heights[i]) if heights[i] is not None else None,651 'wave_direction': float(directions[i]) if directions is not None and i < len(directions) else None,652 'wave_period': float(periods[i]) if periods is not None and i < len(periods) else None,653 'u_velocity': None,654 'v_velocity': None655 }656 regional_points.append(point)657 658 combined_sample_points.extend(regional_points)659 660 # Update bounds for Arctic661 if lats is not None and len(lats) > 0:662 global_lat_min = min(global_lat_min, float(np.min(lats)))663 global_lat_max = max(global_lat_max, float(np.max(lats)))664 global_lon_min = min(global_lon_min, float(np.min(lons)))665 global_lon_max = max(global_lon_max, float(np.max(lons)))666 667 regions_processed.append(region_name)668 logger.info(f"Successfully processed Arctic {region_name}: {len(regional_points)} points")669 670 else:671 # Normal processed data format672 regional_data, _ = result673 if regional_data and 'sample_points' in regional_data:674 # Add regional sample points to global collection675 regional_points = regional_data['sample_points']676 combined_sample_points.extend(regional_points)677 678 # Update global bounds679 grid_info = regional_data.get('grid_info', {})680 if grid_info.get('lat_min') is not None:681 global_lat_min = min(global_lat_min, grid_info['lat_min'])682 global_lat_max = max(global_lat_max, grid_info['lat_max'])683 global_lon_min = min(global_lon_min, grid_info['lon_min'])684 global_lon_max = max(global_lon_max, grid_info['lon_max'])685 686 # Collect wave data for statistics687 for point in regional_points:688 if point.get('wave_height') is not None:689 all_wave_heights.append(point['wave_height'])690 if point.get('wave_direction') is not None:691 all_wave_directions.append(point['wave_direction'])692 if point.get('wave_period') is not None:693 all_wave_periods.append(point['wave_period'])694 695 # Update parameters found696 regional_params = regional_data.get('parameters_found', {})697 if not parameters_found['wave_height']:698 parameters_found['wave_height'] = regional_params.get('wave_height')699 if not parameters_found['wave_direction']:700 parameters_found['wave_direction'] = regional_params.get('wave_direction')701 if not parameters_found['wave_period']:702 parameters_found['wave_period'] = regional_params.get('wave_period')703 if regional_params.get('has_velocity_components'):704 parameters_found['has_velocity_components'] = True705 706 regions_processed.append(region_name)707 logger.info(f"Successfully processed {region_name}: {len(regional_points)} points")708 709 # Clean up temp file710 if os.path.exists(grib_file_path):711 os.unlink(grib_file_path)712 713 except Exception as e:714 logger.error(f"Error processing {region_name} region: {e}")715 # Clean up temp file on error716 if os.path.exists(grib_file_path):717 os.unlink(grib_file_path)718 continue719 720 if not combined_sample_points:721 logger.error("No valid data found in any regional file")722 return None723 724 logger.info(f"Combined data from {len(regions_processed)} regions: {regions_processed}")725 logger.info(f"Total sample points: {len(combined_sample_points)}")726 727 # Create combined global dataset728 combined_data = {729 'timestamp': datetime.utcnow().isoformat(),730 'data_source': f'NOAA_MULTI_REGIONAL_GRIB ({",".join(regions_processed)})',731 'parameters_found': parameters_found,732 'grid_info': {733 'lat_min': float(global_lat_min) if global_lat_min != float('inf') else None,734 'lat_max': float(global_lat_max) if global_lat_max != float('-inf') else None,735 'lon_min': float(global_lon_min) if global_lon_min != float('inf') else None,736 'lon_max': float(global_lon_max) if global_lon_max != float('-inf') else None,737 'regions_included': regions_processed,738 'total_points': len(combined_sample_points)739 },740 'wave_statistics': {741 'max_wave_height': float(max(all_wave_heights)) if all_wave_heights else None,742 'min_wave_height': float(min(all_wave_heights)) if all_wave_heights else None,743 'mean_wave_height': float(np.mean(all_wave_heights)) if all_wave_heights else None,744 'std_wave_height': float(np.std(all_wave_heights)) if all_wave_heights else None745 },746 'sample_points': combined_sample_points747 }748 749 # Add direction/period statistics if available750 if all_wave_directions:751 combined_data['direction_statistics'] = {752 'mean_direction': float(np.mean(all_wave_directions)),753 'direction_std': float(np.std(all_wave_directions))754 }755 756 if all_wave_periods:757 combined_data['period_statistics'] = {758 'max_period': float(max(all_wave_periods)),759 'min_period': float(min(all_wave_periods)),760 'mean_period': float(np.mean(all_wave_periods))761 }762 763 return combined_data764 765 except Exception as e:766 logger.error(f"Error processing multiple regional files: {e}")767 # Clean up any remaining temp files768 for grib_file_path, region_name, _, _ in regional_files:769 if os.path.exists(grib_file_path):770 os.unlink(grib_file_path)771 return None772 773 def _extract_sample_points_with_vectors(self, lats, lons, wave_heights, wave_directions=None, wave_periods=None, num_samples=100):774 """Extract sample points with velocity vector components for visualization"""775 try:776 # Create meshgrid for coordinates777 lon_grid, lat_grid = np.meshgrid(lons, lats)778 779 # Flatten arrays780 flat_lats = lat_grid.flatten()781 flat_lons = lon_grid.flatten()782 flat_waves = wave_heights.flatten()783 784 flat_dirs = None785 flat_periods = None786 787 if wave_directions is not None:788 flat_dirs = wave_directions.flatten()789 if wave_periods is not None:790 flat_periods = wave_periods.flatten()791 792 # Remove NaN values793 valid_mask = ~np.isnan(flat_waves)794 if flat_dirs is not None:795 valid_mask = valid_mask & ~np.isnan(flat_dirs)796 797 valid_lats = flat_lats[valid_mask]798 valid_lons = flat_lons[valid_mask] 799 valid_waves = flat_waves[valid_mask]800 801 if flat_dirs is not None:802 valid_dirs = flat_dirs[valid_mask]803 else:804 valid_dirs = None805 806 if flat_periods is not None:807 valid_periods = flat_periods[valid_mask]808 else:809 valid_periods = None810 811 if len(valid_waves) == 0:812 return []813 814 # Sample points for visualization (to avoid too many points)815 sample_size = min(num_samples, len(valid_waves))816 sample_indices = np.random.choice(len(valid_waves), size=sample_size, replace=False)817 818 sample_points = []819 for idx in sample_indices:820 point = {821 'lat': float(valid_lats[idx]),822 'lon': float(valid_lons[idx]),823 'wave_height': float(valid_waves[idx])824 }825 826 if valid_dirs is not None:827 direction = float(valid_dirs[idx])828 point['wave_direction'] = direction829 830 # Calculate velocity components (u, v) from direction831 # Wave direction is "coming from" in meteorological convention832 # Convert to radians and calculate u,v components833 dir_rad = np.radians(direction)834 # Use wave height as a proxy for wave energy/velocity magnitude835 magnitude = point['wave_height'] * 0.1 # Scale factor for visualization836 837 # Components: u = eastward, v = northward838 # Direction 0° = from North, 90° = from East, etc.839 point['u_component'] = magnitude * np.sin(dir_rad) # eastward component840 point['v_component'] = -magnitude * np.cos(dir_rad) # northward component (negative because "from")841 842 if valid_periods is not None:843 point['wave_period'] = float(valid_periods[idx])844 845 sample_points.append(point)846 847 return sample_points848 849 except Exception as e:850 logger.error(f"Error extracting sample points with vectors: {e}")851 return []852 853 def _extract_sample_points(self, lats, lons, wave_heights, num_samples=100):854 """Legacy method - extract sample points for visualization without vectors"""855 return self._extract_sample_points_with_vectors(lats, lons, wave_heights, None, None, num_samples)856 857 def fetch_global_wave_data(self, forecast_hour=0):858 """Main method to fetch global wave data"""859 try:860 # ECMWF open data doesn't include wave parameters, so use NOAA primarily861 logger.info("Fetching wave data from NOAA WW3 model (ECMWF doesn't provide wave data)...")862 grib_file = None863 model_run = None864 865 # Try NOAA first for wave data866 result = self.fetch_noaa_wave_grib(forecast_hour)867 if result and isinstance(result, list):868 # Multiple regional files downloaded869 regional_files = result870 model_run = regional_files[0][2] if regional_files else None871 grib_file = None # Will process multiple files872 elif result and len(result) == 3:873 # Single file (legacy format)874 grib_file, model_run, actual_forecast_hour = result875 regional_files = None876 else:877 # NOAA failed, try ECMWF as last resort (though it likely won't have wave data)878 logger.info("NOAA failed, trying ECMWF as fallback (unlikely to have wave data)...")879 grib_file = self.fetch_ecmwf_wave_grib(forecast_hour)880 regional_files = None881 model_run = None882 883 if not grib_file and not regional_files:884 logger.error("Both ECMWF and NOAA failed - no real wave data available")885 return None886 887 # Process GRIB file(s)888 if regional_files:889 # Process multiple regional files and combine890 processed_data = self.process_multiple_regional_files(regional_files)891 grib_path = "multiple_regional_files"892 else:893 # Process single GRIB file894 processed_data, grib_path = self.process_grib_file(grib_file)895 896 if processed_data:897 # Add forecast metadata898 processed_data['forecast_info'] = {899 'forecast_hour': forecast_hour,900 'model_run': model_run,901 'forecast_valid_time': (datetime.utcnow() + timedelta(hours=forecast_hour)).isoformat(),902 'is_current': forecast_hour == 0903 }904 905 # Clean up temporary files906 if grib_file and os.path.exists(grib_file):907 os.unlink(grib_file)908 elif regional_files:909 # Files already cleaned up in process_multiple_regional_files910 pass911 912 return processed_data913 914 except Exception as e:915 logger.error(f"Error in fetch_global_wave_data: {e}")916 return None917 918 def fetch_multiple_forecasts(self, forecast_hours=[0, 6, 12, 24, 48]):919 """Fetch multiple forecast time steps"""920 forecasts = {}921 922 for hour in forecast_hours:923 try:924 logger.info(f"Fetching forecast for +{hour} hours...")925 data = self.fetch_global_wave_data(hour)926 if data:927 forecasts[f"f{hour:03d}"] = data928 logger.info(f"Successfully fetched +{hour}h forecast")929 else:930 logger.warning(f"Failed to fetch +{hour}h forecast")931 except Exception as e:932 logger.error(f"Error fetching +{hour}h forecast: {e}")933 continue934 935 return forecasts936 937 def _generate_mock_global_data(self, forecast_hour=0):938 """Generate mock global wave data for testing"""939 logger.info(f"Generating mock global wave data for +{forecast_hour}h forecast...")940 941 # Create a grid of sample points around the world with mock vectors942 sample_points = []943 for lat in range(-60, 61, 20): # Every 20 degrees latitude944 for lon in range(-180, 181, 30): # Every 30 degrees longitude945 # Simulate higher waves in storm-prone areas946 base_height = np.random.uniform(0.5, 2.0)947 if abs(lat) > 40: # Higher latitudes tend to have bigger waves948 base_height += np.random.uniform(0.5, 1.5)949 950 # Generate mock wave direction (random but realistic patterns)951 wave_dir = np.random.uniform(0, 360)952 953 # Calculate velocity components954 magnitude = base_height * 0.1955 dir_rad = np.radians(wave_dir)956 u_comp = magnitude * np.sin(dir_rad)957 v_comp = -magnitude * np.cos(dir_rad)958 959 sample_points.append({960 'lat': float(lat),961 'lon': float(lon),962 'wave_height': round(float(base_height), 2),963 'wave_direction': round(float(wave_dir), 1),964 'wave_period': round(np.random.uniform(4.0, 12.0), 1),965 'u_component': round(float(u_comp), 3),966 'v_component': round(float(v_comp), 3)967 })968 969 return {970 'timestamp': datetime.utcnow().isoformat(),971 'data_source': 'MOCK_GLOBAL_DATA',972 'grid_info': {973 'lat_min': -60.0,974 'lat_max': 60.0,975 'lon_min': -180.0,976 'lon_max': 180.0,977 'lat_resolution': 20.0,978 'lon_resolution': 30.0,979 'grid_shape': [7, 13] # 7 lats x 13 lons980 },981 'wave_statistics': {982 'max_wave_height': max(p['wave_height'] for p in sample_points),983 'min_wave_height': min(p['wave_height'] for p in sample_points),984 'mean_wave_height': np.mean([p['wave_height'] for p in sample_points]),985 'std_wave_height': np.std([p['wave_height'] for p in sample_points])986 },987 'forecast_info': {988 'forecast_hour': forecast_hour,989 'model_run': 'MOCK',990 'forecast_valid_time': (datetime.utcnow() + timedelta(hours=forecast_hour)).isoformat(),991 'is_current': forecast_hour == 0992 },993 'sample_points': sample_points994 }