# import library yang diperlukan
import pandas as pd
import numpy as np
from datetime import *
from netCDF4 import Dataset
from fbprophet import Prophet
import warnings
warnings.filterwarnings('ignore')
File '2010.nc' berisi data rata-rata suhu di beberapa titik koordinat di Asia pada sepanjang tahun 2010.
Akan diprediksi suhu rata-rata di Jakarta pada 30 hari ke depan setelah akhir tahun 2010 (1/1/2011 s/d 30/1/2011) menggunakan FBProphet. Posisi Jakarta adalah di lat: -6.125, long: 106.875
# membaca isi file
data = Dataset('2010.nc')
# menampilkan isi variabel dari dataset
data.variables.keys()
dict_keys(['lon', 'lat', 'time', 'tave', 'rstn'])
# ambil semua data latitude
lat = data.variables['lat'][:]
print(lat[:5])
[-14.875 -14.625 -14.375 -14.125 -13.875]
# mengetahui index posisi latitude=-6.125 (Jakarta) di data latitude
np.where(lat == -6.125)
# index posisi lat ini nanti akan digunakan untuk mendapatkan data suhu di Jakarta
(array([35], dtype=int64),)
ditemukan index posisi latitude Jakarta ada di index ke-35
# ambil semua data longitude
long = data.variables['lon'][:]
print(long[:5])
[60.125 60.375 60.625 60.875 61.125]
# mengetahui index posisi longitude 106.875 (Jakarta) di dalam data longitude
np.where(long == 106.875)
(array([187], dtype=int64),)
diperoleh index posisi longitude Jakarta ada di index ke-187
# ambil semua data time
time = data.variables['time'][:]
print(time[:5])
[ 0. 1440. 2880. 4320. 5760.]
data time di atas dalam satuan menit, dihitung mulai dari tanggal 1/1/2010. Sehingga untuk memudahkan analisis perlu dikonversi ke format date terlebih dahulu.
Misal:
0 menit -> tanggal 1/1/2010
1440 menit -> tanggal 2/1/2010 (1440 menit atau 24 jam atau 1 hari setelah tanggal 1/1/2010)
2880 menit -> tanggal 3/1/2010 (2880 menit atau 48 jam atau 2 hari setelah tanggal 1/1/2010)
dan seterusnya..
# membuat list tanggal dari data time, terhitung mulai tanggal 1/1/2010
dataTime = list(map(lambda x: str(date(2010,1,1) + timedelta(hours=x/60)), time))
# tampilan list tanggal (hanya 5 yang pertama sj yang ditampilkan)
dataTime[:5]
['2010-01-01', '2010-01-02', '2010-01-03', '2010-01-04', '2010-01-05']
# ambil semua data rata-rata suhu
tave = data.variables['tave'][:]
# ambil data rata-rata suhu di Jakarta berdasarkan index posisi latitude dan longitude
temp = tave[:, 35, 187]
# membuat dictionary berisi data time dan data rata-rata suhu di Jakarta
myData = {'time': dataTime, 'temp': temp}
# mengkonversi dictionary ke dataframe
data = pd.DataFrame(myData, columns = ['time','temp'])
# FBProphet mewajibkan penamaan kolom waktunya dengan 'ds', dan kolom nilainya dengan 'y'
# mengubah nama kolom 'time' ke 'ds' dan 'temp' ke 'y'
data = data.rename(columns={'time': 'ds', 'temp': 'y'})
data.head()
ds | y | |
---|---|---|
0 | 2010-01-01 | 28.567158 |
1 | 2010-01-02 | 28.054541 |
2 | 2010-01-03 | 28.693995 |
3 | 2010-01-04 | 29.311541 |
4 | 2010-01-05 | 28.096054 |
# mengeset kolom 'ds' dari dataset sebagai index lalu memplotnya
ax = data.set_index('ds').plot(figsize=(12, 8))
# setting konfigurasi fbprophet pada level konfidensi 95%
# menggunakan analisis seasonality harian, mingguan, bulanan, tahunan
my_model = Prophet(interval_width=0.95, daily_seasonality=True, weekly_seasonality=True, yearly_seasonality=True)
my_model.add_seasonality(name='monthly', period=30.5, fourier_order=3)
<fbprophet.forecaster.Prophet at 0x15b61505c18>
# membuat model
my_model.fit(data)
<fbprophet.forecaster.Prophet at 0x15b61505c18>
# membuat deret tanggal tambahan untuk 30 hari ke depan setelah akhir dataset
future_dates = my_model.make_future_dataframe(periods=30)
future_dates
# jumlah baris data: 365 + 30 = 395
ds | |
---|---|
0 | 2010-01-01 |
1 | 2010-01-02 |
2 | 2010-01-03 |
3 | 2010-01-04 |
4 | 2010-01-05 |
... | ... |
390 | 2011-01-26 |
391 | 2011-01-27 |
392 | 2011-01-28 |
393 | 2011-01-29 |
394 | 2011-01-30 |
395 rows × 1 columns
# memprediksi suhu rata-rata Jakarta 30 hari ke depan (1/1/2011 s/d 30/1/2011)
# yhat_lower dan yhat_upper adalah batas bawah dan atas nilai prediksi
forecast = my_model.predict(future_dates)
forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail(30)
ds | yhat | yhat_lower | yhat_upper | |
---|---|---|---|---|
365 | 2011-01-01 | 27.431444 | 26.057645 | 28.877121 |
366 | 2011-01-02 | 27.165753 | 25.747626 | 28.592039 |
367 | 2011-01-03 | 26.881194 | 25.628534 | 28.259014 |
368 | 2011-01-04 | 26.281654 | 24.862487 | 27.720666 |
369 | 2011-01-05 | 25.954176 | 24.509413 | 27.321670 |
370 | 2011-01-06 | 25.756212 | 24.421035 | 27.256389 |
371 | 2011-01-07 | 25.639396 | 24.311767 | 26.995293 |
372 | 2011-01-08 | 25.640659 | 24.240793 | 27.087974 |
373 | 2011-01-09 | 25.640123 | 24.290092 | 27.119363 |
374 | 2011-01-10 | 25.683511 | 24.284070 | 27.245856 |
375 | 2011-01-11 | 25.365090 | 23.988652 | 26.697579 |
376 | 2011-01-12 | 25.191034 | 23.784440 | 26.545326 |
377 | 2011-01-13 | 24.994737 | 23.605454 | 26.430010 |
378 | 2011-01-14 | 24.767305 | 23.392274 | 26.176372 |
379 | 2011-01-15 | 24.631859 | 23.191522 | 26.032615 |
380 | 2011-01-16 | 24.567544 | 23.207357 | 25.951980 |
381 | 2011-01-17 | 24.693738 | 23.286608 | 26.106070 |
382 | 2011-01-18 | 24.624678 | 23.220019 | 26.034957 |
383 | 2011-01-19 | 24.824873 | 23.502218 | 26.249419 |
384 | 2011-01-20 | 25.039644 | 23.548026 | 26.519674 |
385 | 2011-01-21 | 25.158228 | 23.702214 | 26.690923 |
386 | 2011-01-22 | 25.225765 | 23.761272 | 26.624580 |
387 | 2011-01-23 | 25.195937 | 23.740660 | 26.520023 |
388 | 2011-01-24 | 25.223980 | 23.728253 | 26.502968 |
389 | 2011-01-25 | 25.006968 | 23.596507 | 26.449713 |
390 | 2011-01-26 | 25.107054 | 23.670714 | 26.582997 |
391 | 2011-01-27 | 25.344057 | 23.897160 | 26.698313 |
392 | 2011-01-28 | 25.629057 | 24.201980 | 27.053421 |
393 | 2011-01-29 | 25.966531 | 24.580828 | 27.378171 |
394 | 2011-01-30 | 26.220684 | 24.867594 | 27.565651 |
# memplot titik data asli, dan hasil prediksi
my_model.plot(forecast, uncertainty=True);
Keterangan:
# memplot trend data dari keseluruhan data, serta seasonality: harian, mingguan, bulanan, dan tahunan
my_model.plot_components(forecast);