{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":3354,"databundleVersionId":44362,"sourceType":"competition"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-11-16T22:49:31.301306Z","iopub.execute_input":"2024-11-16T22:49:31.301790Z","iopub.status.idle":"2024-11-16T22:49:31.812408Z","shell.execute_reply.started":"2024-11-16T22:49:31.301742Z","shell.execute_reply":"2024-11-16T22:49:31.810973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install scikit-learn==1.4.2 --quiet","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:31.815538Z","iopub.execute_input":"2024-11-16T22:49:31.816164Z","iopub.status.idle":"2024-11-16T22:49:47.031614Z","shell.execute_reply.started":"2024-11-16T22:49:31.816113Z","shell.execute_reply":"2024-11-16T22:49:47.029959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!unzip /kaggle/input/ams-2014-solar-energy-prediction-contest/gefs_train.zip","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:47.033918Z","iopub.execute_input":"2024-11-16T22:49:47.034498Z","iopub.status.idle":"2024-11-16T22:49:47.042097Z","shell.execute_reply.started":"2024-11-16T22:49:47.034432Z","shell.execute_reply":"2024-11-16T22:49:47.040413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import xarray as xr\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\n\ntqdm.pandas()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:47.045236Z","iopub.execute_input":"2024-11-16T22:49:47.045760Z","iopub.status.idle":"2024-11-16T22:49:47.217411Z","shell.execute_reply.started":"2024-11-16T22:49:47.045691Z","shell.execute_reply":"2024-11-16T22:49:47.215934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load GEFS data","metadata":{}},{"cell_type":"code","source":"# Load the NetCDF file\nfile_path = 'train/dswrf_sfc_latlon_subset_19940101_20071231.nc'\ndata = xr.open_dataset(file_path)\n\n# Display the dataset information\ndata.head()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:47.218967Z","iopub.execute_input":"2024-11-16T22:49:47.219790Z","iopub.status.idle":"2024-11-16T22:49:47.888472Z","shell.execute_reply.started":"2024-11-16T22:49:47.219742Z","shell.execute_reply":"2024-11-16T22:49:47.886992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Select the data variable\ndswrf = data['Downward_Short-Wave_Rad_Flux']\n\nfilled_dswrf = dswrf.interpolate_na(dim='time', method='quadratic', fill_value='extrapolate')\n\n# Compute the average across 'fhour' and 'ens'\naverage_dswrf = filled_dswrf.mean(dim=['fhour', 'ens'])\n\n# Print the result\naverage_dswrf","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:47.890255Z","iopub.execute_input":"2024-11-16T22:49:47.891102Z","iopub.status.idle":"2024-11-16T22:49:50.082720Z","shell.execute_reply.started":"2024-11-16T22:49:47.891052Z","shell.execute_reply":"2024-11-16T22:49:50.081231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stations_df = pd.read_csv(\"/kaggle/input/ams-2014-solar-energy-prediction-contest/station_info.csv\")\nstations_df[\"elon\"] = 360 + stations_df[\"elon\"]\nstations_df.head()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:50.084100Z","iopub.execute_input":"2024-11-16T22:49:50.084661Z","iopub.status.idle":"2024-11-16T22:49:50.217347Z","shell.execute_reply.started":"2024-11-16T22:49:50.084575Z","shell.execute_reply":"2024-11-16T22:49:50.216006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 6))\nplt.imshow(average_dswrf[0],\n           extent=[average_dswrf.lon.min(), average_dswrf.lon.max(),\n                   average_dswrf.lat.min(), average_dswrf.lat.max()],\n           origin='lower', cmap='viridis')\n\nplt.colorbar(label='Average DSWRF')\nplt.scatter(stations_df['elon'], stations_df['nlat'], color='red', label='Stations')\n\nplt.xlabel('Longitude')\nplt.ylabel('Latitude')\nplt.title('Downward Short-Wave Radiation Flux')\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:50.219028Z","iopub.execute_input":"2024-11-16T22:49:50.219412Z","iopub.status.idle":"2024-11-16T22:49:50.780319Z","shell.execute_reply.started":"2024-11-16T22:49:50.219371Z","shell.execute_reply":"2024-11-16T22:49:50.778963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import xarray as xr\nimport pandas as pd\n\ndef get_stations_quantity(quantity_path, data_variable, station_df):\n    data = xr.open_dataset(quantity_path)\n    data_var = data[data_variable]\n    filled_data_var = data_var.interpolate_na(dim='time', method='linear', fill_value='extrapolate')\n    average_data_var = filled_data_var.mean(dim=[dim for dim in ['fhour', 'ens'] if dim in data_var.dims])\n\n    result_df = pd.DataFrame()\n\n    for idx, row in station_df.iterrows():\n        lat = row['nlat']\n        lon = row['elon']\n        stid = row['stid']\n        interpolated_data = average_data_var.interp(lat=lat, lon=lon, method = 'quadratic')\n        temp_df = interpolated_data.to_dataframe().reset_index()\n        temp_df['stid'] = stid\n        temp_df['quantity'] = data_variable\n        temp_df.rename(columns={data_variable: 'value'}, inplace=True)\n\n        result_df = pd.concat([result_df, temp_df], ignore_index=True)\n\n    return result_df\n\ndef combine_datasets(station_df, file_variable_pairs):\n    combined_df = pd.DataFrame()\n\n    for quantity_path, data_variable in tqdm(file_variable_pairs):\n        df = get_stations_quantity(quantity_path, data_variable, station_df)\n        combined_df = pd.concat([combined_df, df], ignore_index=True)\n\n    return combined_df","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:50.782254Z","iopub.execute_input":"2024-11-16T22:49:50.782713Z","iopub.status.idle":"2024-11-16T22:49:50.796344Z","shell.execute_reply.started":"2024-11-16T22:49:50.782659Z","shell.execute_reply":"2024-11-16T22:49:50.794973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define the file and variable pairs\nfile_variable_pairs = [\n    ('train/apcp_sfc_latlon_subset_19940101_20071231.nc', 'Total_precipitation'),\n    ('train/dlwrf_sfc_latlon_subset_19940101_20071231.nc', 'Downward_Long-Wave_Rad_Flux'),\n    ('train/dswrf_sfc_latlon_subset_19940101_20071231.nc', 'Downward_Short-Wave_Rad_Flux'),\n    ('train/pres_msl_latlon_subset_19940101_20071231.nc', 'Pressure'),\n    ('train/pwat_eatm_latlon_subset_19940101_20071231.nc', 'Precipitable_water'),\n    ('train/spfh_2m_latlon_subset_19940101_20071231.nc', 'Specific_humidity_height_above_ground'),\n    ('train/tcdc_eatm_latlon_subset_19940101_20071231.nc', 'Total_cloud_cover'),\n    ('train/tcolc_eatm_latlon_subset_19940101_20071231.nc', 'Total_Column-Integrated_Condensate'),\n    ('train/tmax_2m_latlon_subset_19940101_20071231.nc', 'Maximum_temperature'),\n    ('train/tmin_2m_latlon_subset_19940101_20071231.nc', 'Minimum_temperature'),\n    ('train/tmp_2m_latlon_subset_19940101_20071231.nc', 'Temperature_height_above_ground'),\n    ('train/tmp_sfc_latlon_subset_19940101_20071231.nc', 'Temperature_surface'),\n    ('train/ulwrf_sfc_latlon_subset_19940101_20071231.nc', 'Upward_Long-Wave_Rad_Flux_surface'),\n    ('train/ulwrf_tatm_latlon_subset_19940101_20071231.nc', 'Upward_Long-Wave_Rad_Flux'),\n    ('train/uswrf_sfc_latlon_subset_19940101_20071231.nc', 'Upward_Short-Wave_Rad_Flux'),\n]\n\n# Call the function\ndf_combined = combine_datasets(stations_df, file_variable_pairs)","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:49:50.801299Z","iopub.execute_input":"2024-11-16T22:49:50.801770Z","iopub.status.idle":"2024-11-16T22:52:36.403219Z","shell.execute_reply.started":"2024-11-16T22:49:50.801709Z","shell.execute_reply":"2024-11-16T22:52:36.401763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = df_combined.pivot_table(index=['stid', 'time'], columns=['quantity'], values='value')","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:52:36.405162Z","iopub.execute_input":"2024-11-16T22:52:36.406161Z","iopub.status.idle":"2024-11-16T22:52:43.214271Z","shell.execute_reply.started":"2024-11-16T22:52:36.406098Z","shell.execute_reply":"2024-11-16T22:52:43.212851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train.head()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:52:43.216315Z","iopub.execute_input":"2024-11-16T22:52:43.216948Z","iopub.status.idle":"2024-11-16T22:52:43.247623Z","shell.execute_reply.started":"2024-11-16T22:52:43.216896Z","shell.execute_reply":"2024-11-16T22:52:43.245986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_train = pd.read_csv(\"/kaggle/input/ams-2014-solar-energy-prediction-contest/train.csv\")\ny_train[\"Date\"] = y_train[\"Date\"].apply(lambda x: pd.to_datetime(x, format='%Y%m%d'))","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:52:43.249265Z","iopub.execute_input":"2024-11-16T22:52:43.249770Z","iopub.status.idle":"2024-11-16T22:52:44.139457Z","shell.execute_reply.started":"2024-11-16T22:52:43.249712Z","shell.execute_reply":"2024-11-16T22:52:44.138263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Select a station","metadata":{}},{"cell_type":"code","source":"stat = \"ADAX\"\n\nX_train_stat = X_train[X_train.index.get_level_values('stid') == stat].reset_index()\ny_train_stat = y_train[[\"Date\", stat]]","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:52:44.141072Z","iopub.execute_input":"2024-11-16T22:52:44.141436Z","iopub.status.idle":"2024-11-16T22:52:44.257106Z","shell.execute_reply.started":"2024-11-16T22:52:44.141398Z","shell.execute_reply":"2024-11-16T22:52:44.255566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = X_train_stat.merge(y_train_stat, left_on = \"time\", right_on = \"Date\")","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:52:44.258877Z","iopub.execute_input":"2024-11-16T22:52:44.259411Z","iopub.status.idle":"2024-11-16T22:52:44.269999Z","shell.execute_reply.started":"2024-11-16T22:52:44.259362Z","shell.execute_reply":"2024-11-16T22:52:44.268869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nfrom sklearn.model_selection import train_test_split, GridSearchCV\nfrom sklearn.ensemble import GradientBoostingRegressor\nfrom sklearn.metrics import root_mean_squared_error\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\n\n# Define features and target\nX = df.drop(columns=['stid', 'Date', stat])\ny = df[stat]\n\n# Split the data into training and test sets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n\nX_test_w_time = X_test.copy()\nX_train.drop(columns=[\"time\"], inplace=True)\nX_test.drop(columns=[\"time\"], inplace=True)\n\n# Create a pipeline with scaling and model\npipeline = Pipeline([\n    ('scaler', StandardScaler()),\n    ('model', GradientBoostingRegressor(random_state=42))\n])\n\n# Define parameter grid for grid search\nparam_grid = {\n    'model__n_estimators': [50, 100],\n    'model__max_depth': [3, 5],\n    'model__learning_rate': [0.01, 0.1]\n}\n\n# Set up GridSearchCV\ngrid_search = GridSearchCV(pipeline, param_grid, cv=5, scoring='neg_mean_squared_error', n_jobs=-1)\n\n# Train the model with grid search\ngrid_search.fit(X_train, y_train)\n\n# Best model and parameters\nbest_model = grid_search.best_estimator_\nprint(\"Best params:\", grid_search.best_params_)\n\n# Make predictions\ny_pred = best_model.predict(X_test)\n\n# Calculate RMSE\nrmse = root_mean_squared_error(y_test, y_pred)\nprint(f\"RMSE: {rmse}\")\n\n# Feature importance\nfeature_importances = best_model.named_steps['model'].feature_importances_\nfeatures = X_train.columns\nimportance_df = pd.DataFrame({'Feature': features, 'Importance': feature_importances}).sort_values(by='Importance', ascending=False)\nprint(importance_df)","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:52:44.271575Z","iopub.execute_input":"2024-11-16T22:52:44.272099Z","iopub.status.idle":"2024-11-16T22:53:25.466279Z","shell.execute_reply.started":"2024-11-16T22:52:44.272051Z","shell.execute_reply":"2024-11-16T22:53:25.464700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create a DataFrame with time, ground truth, and predictions\nresult_df = pd.DataFrame({'time': X_test_w_time['time'], 'Ground Truth': y_test, 'Predictions': y_pred})\nresult_df.sort_values('time', inplace = True)\n\n# Plot with the time dimension\nplt.figure(figsize=(14, 7))\nplt.plot(result_df['time'], result_df['Ground Truth'], label='Ground Truth', alpha=0.7)\nplt.plot(result_df['time'], result_df['Predictions'], label='Predictions', alpha=0.3, marker = 'o')\nplt.xlabel('Time')\nplt.ylabel(stat)\nplt.title('Predictions vs Ground Truth over Time')\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-16T22:53:25.468507Z","iopub.execute_input":"2024-11-16T22:53:25.469897Z","iopub.status.idle":"2024-11-16T22:53:26.098487Z","shell.execute_reply.started":"2024-11-16T22:53:25.469841Z","shell.execute_reply":"2024-11-16T22:53:26.096670Z"},"trusted":true},"execution_count":null,"outputs":[]}]}