Feature requested by Issue #26 with additions (#103)

* bug in forecast.run_forecast: 'nwp_source' was not respected when loading the nwp data * adjustment of the test for run_forecast. Test for ICON and GFS NWP. * change data.get_nwp() in a way that the nwp data is available for dates more than 90 days in the past. The hisorical api from open meteo is used for this purpose. As this does not contain any visibility data, the value is set to the maximum of 24000 km. In addition, the 'openmeteo-requests' pip package is used. this is clearer and uses a binary data format instead of json so that the transmission is faster. And so the string tinkering for the api request can be omitted. * remove the request for user input. it is not necessary and makes testing the function unnecessarily complicated. * fix visibility bug and return to a 7 day nwp request just for better comparison * vis as float * cast hole df to float32 since the manual vis add with 24000 was casted as float64 * cast after setting the time index and cast to float64 instead float64 since the old modell also hold dtype float64 --------- Co-authored-by: Jakob Elias Gebler <[email protected]>
openclimatefix · Apr 15, 2024 · 07f8aa2 · 07f8aa2
1 parent 61e6dcd
commit 07f8aa2
Show file tree

Hide file tree

Showing 2 changed files with 89 additions and 60 deletions.
diff --git a/quartz_solar_forecast/data.py b/quartz_solar_forecast/data.py
@@ -9,6 +9,10 @@
 import requests
 import xarray as xr
 
+import openmeteo_requests
+import requests_cache
+from retry_requests import retry
+
 from quartz_solar_forecast.pydantic_models import PVSite
 from quartz_solar_forecast.inverters.enphase import get_enphase_data # Added import for get_enphase_data from /inverters/enphase.py
 
@@ -31,79 +35,101 @@ def get_nwp(site: PVSite, ts: datetime, nwp_source: str = "icon") -> xr.Dataset:
     :return: nwp forecast in xarray
     """
 
+    # Setup the Open-Meteo API client with cache and retry on error
+    cache_session = requests_cache.CachedSession('.cache', expire_after = -1)
+    retry_session = retry(cache_session, retries = 5, backoff_factor = 0.2)
+    openmeteo = openmeteo_requests.Client(session = retry_session)
+
+    # Define the variables we want. Visibility is handled separately after the main request
     variables = [
-        "visibility",
-        "windspeed_10m",
-        "temperature_2m",
-        "precipitation",
-        "shortwave_radiation",
-        "direct_radiation",
-        "cloudcover_low",
-        "cloudcover_mid",
-        "cloudcover_high",
+        "temperature_2m", 
+        "precipitation", 
+        "cloud_cover_low", 
+        "cloud_cover_mid", 
+        "cloud_cover_high", 
+        "wind_speed_10m", 
+        "shortwave_radiation", 
+        "direct_radiation"
     ]
 
     start = ts.date()
     end = start + pd.Timedelta(days=7)
 
-    # Getting NWP, from OPEN METEO
-    url_nwp_source = None
-    if nwp_source == "icon":
-        url_nwp_source = "dwd-icon"
-    elif nwp_source == "gfs":
-        url_nwp_source = "gfs"
+    url = ""
+
+    # check whether the time stamp is more than 3 months in the past
+    if (datetime.now() - ts).days > 90:
+        print("Warning: The requested timestamp is more than 3 months in the past. The weather data are provided by a reanalyse model and not ICON or GFS.")
+
+        # load data from open-meteo Historical Weather API
+        url = "https://archive-api.open-meteo.com/v1/archive"
+
     else:
-        raise Exception(f'Source ({nwp_source}) must be either "icon" or "gfs"')
-
-    # Pull data from the nwp_source provided 
-    url = (
-        f"https://api.open-meteo.com/v1/{url_nwp_source}?"
-        f"latitude={site.latitude}&longitude={site.longitude}"
-        f"&hourly={','.join(variables)}"
-        f"&start_date={start}&end_date={end}"
-    )
-    r = requests.get(url)
-    d = json.loads(r.text)
-
-    # If the nwp_source is ICON, get visibility data from GFS as its not available for icon on Open Meteo
-    if nwp_source == "icon":
-        url = (
-            f"https://api.open-meteo.com/v1/gfs?"
-            f"latitude={site.latitude}&longitude={site.longitude}"
-            f"&hourly=visibility"
-            f"&start_date={start}&end_date={end}"
-        )
-        r_gfs = requests.get(url)
-        d_gfs = json.loads(r_gfs.text)
-
-        # extract visibility data from gfs reponse
-        gfs_visibility_data = d_gfs["hourly"]["visibility"]
-
-        # add visibility to the icon reponse to make a complete json file 
-        d["hourly"]["visibility"] = gfs_visibility_data
 
-    # convert data into xarray
-    df = pd.DataFrame(d["hourly"])
-    df["time"] = pd.to_datetime(df["time"])
-    df = df.rename(
-        columns={
-            "visibility": "vis",
-            "windspeed_10m": "si10",
-            "temperature_2m": "t",
-            "precipitation": "prate",
-            "shortwave_radiation": "dswrf",
-            "direct_radiation": "dlwrf",
-            "cloudcover_low": "lcc",
-            "cloudcover_mid": "mcc",
-            "cloudcover_high": "hcc",
+        # Getting NWP from open meteo weather forecast API by ICON or GFS model within the last 3 months
+        url_nwp_source = None
+        if nwp_source == "icon":
+            url_nwp_source = "dwd-icon"
+        elif nwp_source == "gfs":
+            url_nwp_source = "gfs"
+        else:
+            raise Exception(f'Source ({nwp_source}) must be either "icon" or "gfs"')
+
+        url = f"https://api.open-meteo.com/v1/{url_nwp_source}"
+
+    params = {
+    	"latitude": site.latitude,
+    	"longitude": site.longitude,
+    	"start_date": f"{start}",
+    	"end_date": f"{end}",
+    	"hourly": variables
+    }
+    response = openmeteo.weather_api(url, params=params)
+    hourly = response[0].Hourly()
+
+    hourly_data = {"time": pd.date_range(
+    	start = pd.to_datetime(hourly.Time(), unit = "s", utc = False),
+    	end = pd.to_datetime(hourly.TimeEnd(), unit = "s", utc = False),
+    	freq = pd.Timedelta(seconds = hourly.Interval()),
+    	inclusive = "left"
+    )}
+
+
+    # variables index as in the variables array of the request
+    hourly_data["t"] = hourly.Variables(0).ValuesAsNumpy()
+    hourly_data["prate"] = hourly.Variables(1).ValuesAsNumpy()
+    hourly_data["lcc"] = hourly.Variables(2).ValuesAsNumpy()
+    hourly_data["mcc"] = hourly.Variables(3).ValuesAsNumpy()
+    hourly_data["hcc"] = hourly.Variables(4).ValuesAsNumpy()
+    hourly_data["si10"] = hourly.Variables(5).ValuesAsNumpy()
+    hourly_data["dswrf"] = hourly.Variables(6).ValuesAsNumpy()
+    hourly_data["dlwrf"] = hourly.Variables(7).ValuesAsNumpy()
+
+    # handle visibility 
+    if (datetime.now() - ts).days <= 90:
+        # load data from open-meteo gfs model
+        params = {
+        	"latitude": site.latitude,
+        	"longitude": site.longitude,
+        	"start_date": f"{start}",
+        	"end_date": f"{end}",
+        	"hourly": "visibility"
         }
-    )
+        data_vis_gfs = openmeteo.weather_api("https://api.open-meteo.com/v1/gfs", params=params)[0].Hourly().Variables(0).ValuesAsNumpy()
+        hourly_data["vis"] = data_vis_gfs
+    else:
+        # set to maximum visibility possible
+        hourly_data["vis"] = 24000.0
+
+    df = pd.DataFrame(data = hourly_data)
     df = df.set_index("time")
+    df = df.astype('float64')
+
+    # convert data into xarray
     data_xr = format_nwp_data(df, nwp_source, site)
 
     return data_xr
 
-
 def format_nwp_data(df: pd.DataFrame, nwp_source:str, site: PVSite):
     data_xr = xr.DataArray(
         data=df.values,

diff --git a/requirements.txt b/requirements.txt
@@ -3,4 +3,7 @@ xarray==2022.12.0
 pv-site-prediction==0.1.17
 pydantic==2.6.2
 huggingface_hub== 0.17.3 # only for evaluation
-python-dotenv==1.0.1
+python-dotenv==1.0.1
+openmeteo-requests==1.2.0
+requests-cache==1.2.0
+retry-requests==2.0.0