nakas/ecmwf_open_data_forcast

1#!/usr/bin/env python32"""3Alternative solution using pygrib for polar stereographic GRIB files4"""5 6import numpy as np7import pandas as pd8 9def extract_with_pygrib(grib_file):10    """11    Extract lat/lon and wave data using pygrib library12    """13    try:14        import pygrib15        16        print(f"Opening GRIB file with pygrib: {grib_file}")17        grbs = pygrib.open(grib_file)18        19        # List all messages20        print("Available messages:")21        for i, grb in enumerate(grbs):22            print(f"Message {i+1}: {grb}")23        24        grbs.seek(0)  # Reset to beginning25        26        # Get first message for coordinates27        grb = grbs[1]  # First message28        29        print(f"Grid type: {grb.gridType}")30        print(f"Projection: {getattr(grb, 'projString', 'N/A')}")31        32        # Get lat/lon coordinates33        try:34            lats, lons = grb.latlons()35            print(f"Successfully extracted coordinates")36            print(f"Lat shape: {lats.shape}, range: {lats.min():.3f} to {lats.max():.3f}")37            print(f"Lon shape: {lons.shape}, range: {lons.min():.3f} to {lons.max():.3f}")38            39            # Get data values40            data = grb.values41            print(f"Data shape: {data.shape}, range: {data.min():.3f} to {data.max():.3f}")42            43            # Create flat arrays for lat, lon, data44            lats_flat = lats.flatten()45            lons_flat = lons.flatten()46            data_flat = data.flatten()47            48            # Remove any invalid data points49            valid_mask = ~np.isnan(data_flat) & ~np.isnan(lats_flat) & ~np.isnan(lons_flat)50            51            result_df = pd.DataFrame({52                'latitude': lats_flat[valid_mask],53                'longitude': lons_flat[valid_mask],54                'wave_data': data_flat[valid_mask]55            })56            57            print(f"Created DataFrame with {len(result_df)} valid points")58            59            grbs.close()60            return result_df61            62        except Exception as e:63            print(f"Error extracting lat/lon with pygrib: {e}")64            65            # Try alternative coordinate extraction66            try:67                # Some GRIB files store coordinates differently68                if hasattr(grb, 'latitudeOfFirstGridPointInDegrees'):69                    lat_first = grb.latitudeOfFirstGridPointInDegrees70                    lon_first = grb.longitudeOfFirstGridPointInDegrees71                    print(f"First grid point: lat={lat_first}, lon={lon_first}")72                73                if hasattr(grb, 'Ni') and hasattr(grb, 'Nj'):74                    ni, nj = grb.Ni, grb.Nj75                    print(f"Grid dimensions: {ni} x {nj}")76                    77                # Get projection parameters78                if hasattr(grb, 'projParams'):79                    print(f"Projection parameters: {grb.projParams}")80                    81            except Exception as e2:82                print(f"Error getting grid parameters: {e2}")83            84            grbs.close()85            return None86            87    except ImportError:88        print("pygrib not available. Install with: pip install pygrib")89        return None90    except Exception as e:91        print(f"Error with pygrib: {e}")92        return None93 94def manual_polar_stereographic_projection(grib_file):95    """96    Manually implement polar stereographic coordinate calculation97    """98    try:99        import pygrib100        101        grbs = pygrib.open(grib_file)102        grb = grbs[1]103        104        # Get polar stereographic parameters105        params = {}106        attr_names = [107            'latitudeOfFirstGridPointInDegrees',108            'longitudeOfFirstGridPointInDegrees', 109            'DxInMetres',110            'DyInMetres',111            'orientationOfTheGridInDegrees',112            'latitudeWhereDxAndDyAreSpecifiedInDegrees',113            'Ni', 'Nj'114        ]115        116        for attr in attr_names:117            if hasattr(grb, attr):118                params[attr] = getattr(grb, attr)119                print(f"{attr}: {params[attr]}")120        121        # Manual coordinate calculation for polar stereographic122        if 'Ni' in params and 'Nj' in params:123            ni, nj = int(params['Ni']), int(params['Nj'])124            125            # Get grid spacing126            dx = params.get('DxInMetres', 25000)  # Default 25km127            dy = params.get('DyInMetres', 25000)128            129            # Create coordinate arrays130            x = np.arange(ni) * dx131            y = np.arange(nj) * dy132            133            X, Y = np.meshgrid(x, y)134            135            # Convert from polar stereographic to lat/lon136            # This is a simplified version - you may need a proper projection library137            lat_origin = params.get('latitudeWhereDxAndDyAreSpecifiedInDegrees', 90.0)138            lon_origin = params.get('orientationOfTheGridInDegrees', 0.0)139            140            print(f"Grid origin: lat={lat_origin}, lon={lon_origin}")141            print(f"Grid spacing: dx={dx}m, dy={dy}m")142            print(f"Grid size: {ni} x {nj}")143            144            # For proper conversion, you'd use pyproj or similar145            try:146                from pyproj import Proj, transform147                148                # Define polar stereographic projection149                proj_polar = Proj(150                    proj='stere',151                    lat_0=lat_origin,152                    lon_0=lon_origin,153                    lat_ts=lat_origin,154                    ellps='sphere'155                )156                157                proj_latlon = Proj(proj='latlong', ellps='sphere')158                159                # Transform coordinates160                lons, lats = transform(proj_polar, proj_latlon, X.flatten(), Y.flatten())161                162                lats = np.array(lats).reshape(nj, ni)163                lons = np.array(lons).reshape(nj, ni)164                165                # Get data166                data = grb.values167                168                # Create DataFrame169                result_df = pd.DataFrame({170                    'latitude': lats.flatten(),171                    'longitude': lons.flatten(),172                    'wave_data': data.flatten()173                })174                175                # Remove invalid points176                valid_mask = ~np.isnan(result_df['wave_data'])177                result_df = result_df[valid_mask]178                179                print(f"Successfully calculated coordinates for {len(result_df)} points")180                181                grbs.close()182                return result_df183                184            except ImportError:185                print("pyproj not available for coordinate transformation")186                grbs.close()187                return None188        189        grbs.close()190        return None191        192    except Exception as e:193        print(f"Error in manual projection: {e}")194        return None195 196# Example usage197if __name__ == "__main__":198    grib_file = "/tmp/tmpr004q4kw.grib2"  # Your file path199    200    print("=== Trying pygrib extraction ===")201    result1 = extract_with_pygrib(grib_file)202    203    if result1 is not None:204        print(f"Success with pygrib! {len(result1)} points extracted")205        print(result1.head())206    else:207        print("=== Trying manual polar stereographic ===")208        result2 = manual_polar_stereographic_projection(grib_file)209        210        if result2 is not None:211            print(f"Success with manual projection! {len(result2)} points extracted")212            print(result2.head())213        else:214            print("All methods failed. Consider using wgrib2 conversion first.")215