Download and aggregate ECMWF Open Data 6 hour forecasts to daily timestep and create map-stacks for the next 10 days for the African continent¶
Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands
The ECMWF Open Data service provides free access to selected forecast products from the ECMWF Integrated Forecasting System (IFS). It is intended for operational forecasting, research, and educational use. For further information see: https://www.ecmwf.int/en/forecasts/datasets/open-data.
The main characteristics are:
| Feature | ECMWF Open Data |
|---|---|
| Data type | Operational forecasts |
| Forecast frequency | 00 UTC and 12 UTC |
| Forecast horizon | Up to 10 days (higher frequency initially, then longer intervals) |
| Spatial resolution | Approximately 0.25° (~28 km) |
| Variables | Temperature, precipitation, wind, pressure, humidity, radiation, etc. |
| Format | GRIB2 |
| Access | HTTP and AWS S3 |
In [1]:
#!pip install ecmwf-opendata #https://pypi.org/project/ecmwf-opendata/
In [2]:
import xarray as xr
import numpy as np
import pandas as pd
import os
from ecmwf.opendata import Client
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import ilwis
import warnings
warnings.filterwarnings("ignore")
In [3]:
#downloaded data and newly created animations in same folder as notebook, change accordingly!
work_dir = os.getcwd()
In [4]:
#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)
d:\jupyter\notebook_scripts\Special\ECMWF_OpenData
In [5]:
#Link to the client - source is either the name of server to contact or a fully qualified URL.
#Possible values are ecmwf to access ECMWF's servers, aws for data hosted by Amazon Web Services,
#or azure to access data hosted on Microsoft's Azure. Default is ecmwf.
#Known sources are ['ecmwf', 'azure', 'aws', 'ecmwf-esuites']
#client = Client(source="aws", model="ifs", resol="0p25") #0.25 degree
client = Client(source="ecmwf", model="ifs", resol="0p25") #0.25 degree
In [6]:
#retrieve the 10 day forecast, 6 hourly interval
steps = list(range(0, 240, 6))
print(steps)
[0, 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84, 90, 96, 102, 108, 114, 120, 126, 132, 138, 144, 150, 156, 162, 168, 174, 180, 186, 192, 198, 204, 210, 216, 222, 228, 234]
In [7]:
# uncomment lines below to retrieve latest forecast - all parameters at once
# client.retrieve(
# param=["tp", "msl", "10u", "10v", "2t", "ro"],
# step=steps,
# type="fc",
# target="latest_fc_ECMWF.grib2"
# )
In [8]:
target="latest_fc_ECMWF.grib2"
print(str(target))
latest_fc_ECMWF.grib2
In [9]:
t2m = xr.open_dataset(target, engine="cfgrib", filter_by_keys={"shortName": "2t"}) # K
msl = xr.open_dataset(target, engine="cfgrib", filter_by_keys={"shortName": "msl"}) # Pa
tp = xr.open_dataset(target, engine="cfgrib", filter_by_keys={"shortName": "tp"}) # accumulated precipitation
ro = xr.open_dataset(target, engine="cfgrib", filter_by_keys={"shortName": "ro"}) # accumulated runoff
u10 = xr.open_dataset(target, engine="cfgrib", filter_by_keys={"shortName": "10u"}) # m/s
v10 = xr.open_dataset(target, engine="cfgrib", filter_by_keys={"shortName": "10v"}) # m/s
In [10]:
#data check of selected parameter
print(t2m)
<xarray.Dataset> Size: 166MB
Dimensions: (step: 40, latitude: 721, longitude: 1440)
Coordinates:
time datetime64[ns] 8B ...
* step (step) timedelta64[ns] 320B 0 days 00:00:00 ... 9 days...
heightAboveGround float64 8B ...
* latitude (latitude) float64 6kB 90.0 89.75 89.5 ... -89.75 -90.0
* longitude (longitude) float64 12kB -180.0 -179.8 ... 179.5 179.8
valid_time (step) datetime64[ns] 320B ...
Data variables:
t2m (step, latitude, longitude) float32 166MB ...
Attributes:
GRIB_edition: 2
GRIB_centre: ecmf
GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts
GRIB_subCentre: 0
Conventions: CF-1.7
institution: European Centre for Medium-Range Weather Forecasts
history: 2026-07-13T14:34 GRIB to CDM+CF via cfgrib-0.9.1...
In [11]:
valid_time = t2m.time + t2m.step
In [12]:
for ds in [t2m, msl, tp, ro, u10, v10]:
ds.coords["valid_time"] = valid_time
In [13]:
print(t2m.valid_time.values)
['2026-07-13T00:00:00.000000000' '2026-07-13T06:00:00.000000000' '2026-07-13T12:00:00.000000000' '2026-07-13T18:00:00.000000000' '2026-07-14T00:00:00.000000000' '2026-07-14T06:00:00.000000000' '2026-07-14T12:00:00.000000000' '2026-07-14T18:00:00.000000000' '2026-07-15T00:00:00.000000000' '2026-07-15T06:00:00.000000000' '2026-07-15T12:00:00.000000000' '2026-07-15T18:00:00.000000000' '2026-07-16T00:00:00.000000000' '2026-07-16T06:00:00.000000000' '2026-07-16T12:00:00.000000000' '2026-07-16T18:00:00.000000000' '2026-07-17T00:00:00.000000000' '2026-07-17T06:00:00.000000000' '2026-07-17T12:00:00.000000000' '2026-07-17T18:00:00.000000000' '2026-07-18T00:00:00.000000000' '2026-07-18T06:00:00.000000000' '2026-07-18T12:00:00.000000000' '2026-07-18T18:00:00.000000000' '2026-07-19T00:00:00.000000000' '2026-07-19T06:00:00.000000000' '2026-07-19T12:00:00.000000000' '2026-07-19T18:00:00.000000000' '2026-07-20T00:00:00.000000000' '2026-07-20T06:00:00.000000000' '2026-07-20T12:00:00.000000000' '2026-07-20T18:00:00.000000000' '2026-07-21T00:00:00.000000000' '2026-07-21T06:00:00.000000000' '2026-07-21T12:00:00.000000000' '2026-07-21T18:00:00.000000000' '2026-07-22T00:00:00.000000000' '2026-07-22T06:00:00.000000000' '2026-07-22T12:00:00.000000000' '2026-07-22T18:00:00.000000000']
In [14]:
#get info on number of timesteps
t2m = t2m.swap_dims({"step": "valid_time"})
msl = msl.swap_dims({"step": "valid_time"})
tp = tp.swap_dims({"step": "valid_time"})
ro = ro.swap_dims({"step": "valid_time"})
u10 = u10.swap_dims({"step": "valid_time"})
v10 = v10.swap_dims({"step": "valid_time"})
In [15]:
print(t2m.dims)
FrozenMappingWarningOnValuesAccess({'valid_time': 40, 'latitude': 721, 'longitude': 1440})
In [16]:
#map extent
lat_min, lat_max = -38, 40
lon_min, lon_max = -26, 52
In [17]:
#create subset
t2m = t2m.sel(latitude=slice(lat_max, lat_min), longitude=slice(lon_min, lon_max))
msl = msl.sel(latitude=slice(lat_max, lat_min), longitude=slice(lon_min, lon_max))
tp = tp.sel(latitude=slice(lat_max, lat_min), longitude=slice(lon_min, lon_max))
ro = ro.sel(latitude=slice(lat_max, lat_min), longitude=slice(lon_min, lon_max))
u10 = u10.sel(latitude=slice(lat_max, lat_min), longitude=slice(lon_min, lon_max))
v10 = v10.sel(latitude=slice(lat_max, lat_min), longitude=slice(lon_min, lon_max))
In [18]:
#convert from Kelvin to Celsius and resample to daily timestep
t2m_daily_mean = t2m["t2m"].resample(valid_time="1D").mean() - 273.15
t2m_daily_max = t2m["t2m"].resample(valid_time="1D").max() - 273.15
t2m_daily_min = t2m["t2m"].resample(valid_time="1D").min() - 273.15
In [19]:
t2m_daily_max
Out[19]:
<xarray.DataArray 't2m' (valid_time: 10, latitude: 313, longitude: 313)> Size: 4MB
array([[[19.900757 , 19.994507 , 19.994507 , ..., 27.775757 ,
27.400757 , 27.182007 ],
[20.182007 , 20.119507 , 20.150757 , ..., 27.650757 ,
27.463257 , 27.057007 ],
[20.213257 , 20.213257 , 20.182007 , ..., 27.557007 ,
27.057007 , 26.932007 ],
...,
[11.025757 , 10.869507 , 10.932007 , ..., 10.838257 ,
10.744507 , 10.838257 ],
[10.775757 , 10.713257 , 10.744507 , ..., 10.932007 ,
10.869507 , 10.775757 ],
[10.713257 , 10.650757 , 10.619507 , ..., 10.963257 ,
10.869507 , 10.463257 ]],
[[19.675018 , 19.612518 , 19.771881 , ..., 27.115631 ,
26.959381 , 27.048187 ],
[19.768768 , 19.800018 , 19.831268 , ..., 27.209381 ,
27.053131 , 26.860687 ],
[19.956268 , 19.956268 , 19.925018 , ..., 27.240631 ,
27.021881 , 26.923187 ],
...
[16.00418 , 16.06668 , 16.06668 , ..., 14.951141 ,
15.138641 , 15.263641 ],
[15.879181 , 15.941681 , 15.816681 , ..., 15.201141 ,
15.138641 , 15.138641 ],
[15.754181 , 15.691681 , 15.566681 , ..., 15.263641 ,
15.263641 , 15.138641 ]],
[[21.616058 , 21.397308 , 21.334808 , ..., 30.137817 ,
29.700317 , 29.637817 ],
[21.459808 , 21.366058 , 21.397308 , ..., 29.491058 ,
29.887817 , 29.887817 ],
[21.522308 , 21.459808 , 21.428558 , ..., 29.334808 ,
29.397308 , 29.366058 ],
...,
[13.169464 , 13.356964 , 13.481964 , ..., 16.084808 ,
16.084808 , 16.116058 ],
[13.044464 , 13.044464 , 12.919464 , ..., 16.209808 ,
16.147308 , 16.116058 ],
[12.731964 , 12.856964 , 12.981964 , ..., 16.116058 ,
16.147308 , 16.084808 ]]], shape=(10, 313, 313), dtype=float32)
Coordinates:
time datetime64[ns] 8B 2026-07-13
heightAboveGround float64 8B 2.0
* latitude (latitude) float64 3kB 40.0 39.75 39.5 ... -37.75 -38.0
* longitude (longitude) float64 3kB -26.0 -25.75 -25.5 ... 51.75 52.0
* valid_time (valid_time) datetime64[ns] 80B 2026-07-13 ... 2026-07-22In [20]:
#convert from Pascal to hPa and resample to (mean) daily timestep
msl_daily = msl["msl"].resample(valid_time="1D").mean() / 100
In [21]:
#process wind speed and directions and resample to (mean) daily timestep
u_daily = u10["u10"].resample(valid_time="1D").mean()
v_daily = v10["v10"].resample(valid_time="1D").mean()
ws_daily = np.sqrt(u_daily**2 + v_daily**2)
wd_daily = (180 + np.degrees(np.arctan2(u_daily, v_daily))) % 360
In [22]:
#retrieve accumulated precipitation
#The .diff() function computes the difference between consecutive time steps along the valid_time dimension
tp_interval = tp["tp"].diff(dim="valid_time")
#Since .diff() removes the first time step, pad() adds it back
#valid_time=(1,0), means add 1 value before the existing data, add 0 values after. And constant_values=0, ensures that new first value is filled with 0
tp_interval = tp_interval.pad(valid_time=(1,0), constant_values=0)
#After padding, the coordinate labels are no longer exactly the same as the original dataset.
#This line replaces the time coordinates with the original ones
tp_interval = tp_interval.assign_coords(valid_time=tp.valid_time.values)
#Print resulting timesteps
print(tp_interval.valid_time.values[:10])
['2026-07-13T00:00:00.000000000' '2026-07-13T06:00:00.000000000' '2026-07-13T12:00:00.000000000' '2026-07-13T18:00:00.000000000' '2026-07-14T00:00:00.000000000' '2026-07-14T06:00:00.000000000' '2026-07-14T12:00:00.000000000' '2026-07-14T18:00:00.000000000' '2026-07-15T00:00:00.000000000' '2026-07-15T06:00:00.000000000']
In [23]:
#convert from meter to mm and aggregate precipitation values (sum) to daily timestep
tp_daily = tp_interval.resample(valid_time="1D").sum() * 1000
tp_daily
Out[23]:
<xarray.DataArray 'tp' (valid_time: 10, latitude: 313, longitude: 313)> Size: 4MB
array([[[0.0000000e+00, 7.6293945e-03, 1.5258789e-02, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
[0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
[0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
...,
[1.6021729e-01, 2.8991699e-01, 1.2207031e-01, ...,
4.8828125e-01, 5.2642822e-01, 4.3487549e-01],
[2.0599365e-01, 1.6021729e-01, 2.3651123e-01, ...,
6.1798096e-01, 7.0190430e-01, 6.2561035e-01],
[2.2888184e-01, 1.8310547e-01, 1.0681152e-01, ...,
6.3323975e-01, 7.0190430e-01, 5.2642822e-01]],
[[2.8228760e-01, 2.8228760e-01, 3.7384033e-01, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
[2.8991699e-01, 2.5939941e-01, 2.9754639e-01, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
[2.8228760e-01, 2.4414062e-01, 2.3651123e-01, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
...
2.5939941e-01, 2.7465820e-01, 1.9836426e-01],
[2.2613525e+01, 2.1881104e+01, 2.5848389e+01, ...,
3.6621094e-01, 3.6621094e-01, 3.2043457e-01],
[2.9800415e+01, 2.9891968e+01, 2.5451660e+01, ...,
3.5095215e-01, 3.8146973e-01, 4.1198730e-01]],
[[2.0599365e+01, 1.8676758e+01, 1.2512207e+01, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
[1.4282227e+01, 1.3732910e+01, 8.7890625e+00, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
[8.2397461e+00, 6.5917969e+00, 4.3334961e+00, ...,
0.0000000e+00, 0.0000000e+00, 0.0000000e+00],
...,
[3.6621094e+00, 2.8991699e+00, 2.8381348e+00, ...,
5.4931641e-01, 6.7138672e-01, 1.5258789e-01],
[2.6245117e+00, 2.7770996e+00, 2.0141602e+00, ...,
7.9345703e-01, 9.1552734e-01, 9.7656250e-01],
[2.3498535e+00, 1.5869141e+00, 1.4343262e+00, ...,
6.4086914e-01, 7.0190430e-01, 8.5449219e-01]]],
shape=(10, 313, 313), dtype=float32)
Coordinates:
time datetime64[ns] 8B 2026-07-13
surface float64 8B ...
* latitude (latitude) float64 3kB 40.0 39.75 39.5 ... -37.75 -38.0
* longitude (longitude) float64 3kB -26.0 -25.75 -25.5 ... 51.75 52.0
heightAboveGround float64 8B 2.0
* valid_time (valid_time) datetime64[ns] 80B 2026-07-13 ... 2026-07-22In [24]:
ro_interval = ro["ro"].diff(dim="valid_time")
ro_interval = ro_interval.pad(valid_time=(1, 0), constant_values=0)
ro_interval = ro_interval.assign_coords(valid_time=ro.valid_time.values)
#convert from meter to mm and aggregate total runoff values to daily timestep
ro_daily = ro_interval.resample(valid_time="1D").sum() * 1000
In [25]:
ro_daily
Out[25]:
<xarray.DataArray 'ro' (valid_time: 10, latitude: 313, longitude: 313)> Size: 4MB
array([[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
...
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]],
[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]],
shape=(10, 313, 313), dtype=float32)
Coordinates:
time datetime64[ns] 8B 2026-07-13
surface float64 8B ...
* latitude (latitude) float64 3kB 40.0 39.75 39.5 ... -37.75 -38.0
* longitude (longitude) float64 3kB -26.0 -25.75 -25.5 ... 51.75 52.0
heightAboveGround float64 8B 2.0
* valid_time (valid_time) datetime64[ns] 80B 2026-07-13 ... 2026-07-22In [26]:
# =========================================================
# KEEP FIRST 10 DAYS
# =========================================================
tp_daily_mm = tp_daily.isel(valid_time=slice(0, 10))
ro_daily_mm = ro_daily.isel(valid_time=slice(0, 10))
t2m_daily_mean_c = t2m_daily_mean.isel(valid_time=slice(0, 10))
t2m_daily_min_c = t2m_daily_min.isel(valid_time=slice(0, 10))
t2m_daily_max_c = t2m_daily_max.isel(valid_time=slice(0, 10))
msl_daily_hpa = msl_daily.isel(valid_time=slice(0, 10))
wind_speed_daily = ws_daily.isel(valid_time=slice(0, 10))
wind_dir_daily = wd_daily.isel(valid_time=slice(0, 10))
In [27]:
# =========================================================
# SUMMARY
# =========================================================
print("Daily precipitation:", tp_daily_mm.shape)
print("Daily runoff:", ro_daily_mm.shape)
print("Daily temperature mean:", t2m_daily_mean_c.shape)
print("Daily temperature min:", t2m_daily_min_c.shape)
print("Daily temperature max:", t2m_daily_max_c.shape)
print("Daily pressure:", msl_daily_hpa.shape)
print("Daily wind speed:", wind_speed_daily.shape)
print("Daily wind direction:", wind_dir_daily.shape)
Daily precipitation: (10, 313, 313) Daily runoff: (10, 313, 313) Daily temperature mean: (10, 313, 313) Daily temperature min: (10, 313, 313) Daily temperature max: (10, 313, 313) Daily pressure: (10, 313, 313) Daily wind speed: (10, 313, 313) Daily wind direction: (10, 313, 313)
In [28]:
# First day
temp_day1 = t2m_daily_mean_c.isel(valid_time=0)
prec_day1 = tp_daily_mm.isel(valid_time=0)
prec_day1_masked = np.ma.masked_less(prec_day1, 0.1) #mask no or very low pcp
# Create figure with two subplots
fig, axes = plt.subplots(
1, 2,
figsize=(14, 6),
subplot_kw={"projection": ccrs.PlateCarree()}
)
# ==========================
# Temperature
# ==========================
ax = axes[0]
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
ax.coastlines(linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
ax.add_feature(cfeature.RIVERS, linewidth=0.5)
ax.add_feature(cfeature.LAKES)
im1 = ax.pcolormesh(
temp_day1.longitude,
temp_day1.latitude,
temp_day1,
cmap="RdYlBu_r",
shading="auto",
transform=ccrs.PlateCarree()
)
cbar1 = plt.colorbar(im1, ax=ax, shrink=0.7, pad=0.03)
cbar1.set_label("Temperature (°C)")
# ==========================
# Precipitation
# ==========================
ax = axes[1]
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
ax.coastlines(linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
ax.add_feature(cfeature.RIVERS, linewidth=0.5)
ax.add_feature(cfeature.LAKES)
im2 = ax.pcolormesh(
prec_day1.longitude,
prec_day1.latitude,
prec_day1_masked,
cmap="jet",
shading="auto",
vmin=0,
vmax=25,
transform=ccrs.PlateCarree()
)
cbar2 = plt.colorbar(im2, ax=ax, shrink=0.7, pad=0.03)
cbar2.set_label("Precipitation (mm/day)")
# Overall title
date = str(temp_day1.valid_time.values)[:10]
fig.suptitle(f"ECMWF Daily Weather Fields\n{date}", fontsize=15)
plt.tight_layout()
plt.show()
In [29]:
grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=-26 40 52 -38 , gridsize=313 313, cornerofcorners=yes')
In [30]:
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-100, 100.0, 0))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(161, 213,10))
rcNew.setGeoReference(grf_LL)
rcNew.setDataDef(dfNum)
print(rcNew.size())
Size(313, 313, 10)
Export the daily forecasts¶
In [31]:
#T-Mean
mapstack_out = np.array([t2m_daily_mean_c]).flatten()
print(mapstack_out.shape)
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('T_mean.tif', 'GTiff', 'gdal')
(979690,) Size(313, 313, 10)
In [32]:
#T-Min
mapstack_out = np.array([t2m_daily_min_c]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('T_min.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
In [33]:
#T-Max
mapstack_out = np.array([t2m_daily_max_c]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('T_max.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
In [34]:
#Daily precipitation
mapstack_out = np.array([tp_daily_mm]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('Prec.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
In [35]:
#Daily runoff
mapstack_out = np.array([ro_daily_mm]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('Runoff.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
In [36]:
#Daily pressure
mapstack_out = np.array([msl_daily_hpa]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('Pressure.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
In [37]:
#Daily wind speed
mapstack_out = np.array([wind_speed_daily]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('W_speed.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
In [38]:
#Daily wind direction
mapstack_out = np.array([wind_dir_daily]).flatten()
rcNew.array2raster(mapstack_out)
print(rcNew.size())
rcNew.store('W_dir.tif', 'GTiff', 'gdal')
Size(313, 313, 10)
Check your results using ILWIS desktop, import the geotif files and display the maplists as animated sequences.
In [ ]: