[1] "Hello, R!"Analisis dan Visualisasi Data Iklim dengan R & Python
Klimatologi Terapan - IPB University
Jenis Data Iklim
- Data iklim memiliki berbagai bentuk, mulai dari data observasi, hasil model (biasanya berbentuk grid), serta data citra satelit
- Format tipe data iklim: spreadsheet/database (seperti Excel, CSV, TSV, SQLite, …), netCDF, HDF, GRIB, dan Zarr.
- Dimensi data iklim: waktu (t), lokasi (x, y, z), dan variabel (pr, t2m, u10, v10, …)
- Semakin berkembangnya teknologi komputasi, data iklim saat ini semakin banyak sehingga ukuran file datanya juga semakin besar. Contoh data iklim tersedia gratis: Global Forecast System (GFS), ECMWF IFS, dan ERA5
Tip
Untuk mengubah format data iklim GRIB menjadi netCDF, Anda dapat menggunakan software Climate Data Operator (CDO) dan perintahnya adalah cdo -f nc copy nama-file.grib nama-file.nc
Spasialisasi Data Iklim (2)
library(gstat)
library(raster)
library(sf)
library(tmap)
tmap_options(check.and.fix = TRUE)
# Filter Bulan: Januari 2019 & 2022
dats_monthly_jan <- dats_monthly %>% filter(month(Bulan) == 1 & year(Bulan) %in% c(2019, 2022))
dats_jan_2019 <- dats_monthly_jan %>% filter(Tahun == 2019)
dats_jan_2022 <- dats_monthly_jan %>% filter(Tahun == 2022)
# Convert dataframe to SpatialPointsDataFrames
coordinates(dats_jan_2019) <- ~lon + lat
coordinates(dats_jan_2022) <- ~lon + lat
# Make an empty grid with the same extent as the data and resolution set to 0.05
res <- 0.05
empty_grd <- sf::st_make_grid(dats_jan_2019_sf, what = 'centers', cellsize = 0.05, crs = 4328)
# Do the IDW calculation
power <- 1
idw_jan2019 <- idw(RR ~ 1, dats_jan_2019_sf, newdata = empty_grd, idp = power)
idw_jan2022 <- idw(RR ~ 1, dats_jan_2022_sf, newdata = empty_grd, idp = power)
# Assign coordinates to the IDW result
idw_jan2019 <- idw_jan2019 %>%
mutate(lon = st_coordinates(idw_jan2019)[, 1],
lat = st_coordinates(idw_jan2019)[, 2])
idw_jan2022 <- idw_jan2022 %>%
mutate(lon = st_coordinates(idw_jan2022)[, 1],
lat = st_coordinates(idw_jan2022)[, 2])
# Sort lat descending and lon ascending
# This because of raster package requires the data to be sorted in this way
# From the up left corner to the bottom right corner
idw_jan2019 <- idw_jan2019 %>%
arrange(desc(lat)) %>%
arrange(lon)
idw_jan2022 <- idw_jan2022 %>%
arrange(desc(lat)) %>%
arrange(lon)
# convert to matrix which has ncell = lon * lat
lon <- unique(idw_jan2019$lon); nlon <- length(lon)
lat <- unique(idw_jan2019$lat); nlat <- length(lat)
idw_mat_jan2019 <- matrix(idw_jan2019$var1.pred, nrow = nlat, ncol = nlon)
idw_mat_jan2022 <- matrix(idw_jan2022$var1.pred, nrow = nlat, ncol = nlon)
# Convert to raster
idw_raster_jan2019 <- raster(idw_mat_jan2019, xmn = min(lon), xmx = max(lon),
ymn = min(lat), ymx = max(lat))
idw_raster_jan2022 <- raster(idw_mat_jan2022, xmn = min(lon), xmx = max(lon),
ymn = min(lat), ymx = max(lat))
# Mask the raster with the boundary of Pulau Jawa
p_jawa <- st_read("data/jawa.geojson")
stck <- stack(idw_raster_jan2019, idw_raster_jan2022)
stck <- mask(stck, p_jawa)
names(stck) <- c("IDW Jan_2019", "IDW Jan_2022")
# Plot the result
idw_plt <- tm_shape(stck) +
tm_raster(style = 'cont') +
tm_shape(p_jawa) +
tm_borders() +
tm_facets(nrow = 2)
idw_plt
import geopandas as gpd
import numpy as np
import xarray as xr
import rioxarray as rxr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LinearSegmentedColormap
# Filter for January 2019 and 2022
dats_monthly_jan = dats_monthly[(dats_monthly['Bulan'].dt.month == 1) & (dats_monthly['Bulan'].dt.year.isin([2019, 2022]))]
dats_monthly_jan_2019 = dats_monthly_jan[dats_monthly_jan['Bulan'].dt.year == 2019]
dats_monthly_jan_2022 = dats_monthly_jan[dats_monthly_jan['Bulan'].dt.year == 2022]
# Define grid space
grid_space = 0.05
# Define grid_lon and grid_lat based on the min and max of lon and lat
grid_lon = np.arange(dats_monthly_jan['lon'].min(), dats_monthly_jan['lon'].max(), grid_space)
grid_lat = np.arange(dats_monthly_jan['lat'].min(), dats_monthly_jan['lat'].max(), grid_space)
# Combine grid_lon and grid_lat into a meshgrid
all_lats = np.meshgrid(grid_lon, grid_lat)[1].ravel()
all_lons = np.meshgrid(grid_lon, grid_lat)[0].ravel()
# Make empty data frame and assign all_lats and all_lons to the data frame
itrp = pd.DataFrame()
itrp['lat'] = all_lats
itrp['lon'] = all_lons
# Do the IDW with sklearn.neighbors.KNeighborsRegressor
from sklearn.neighbors import KNeighborsRegressor
# Define the number of neighbors
n_neighbors = 10
# Define the power parameter
power = 1
# Define the KNeighborsRegressor
knn = KNeighborsRegressor(n_neighbors=n_neighbors, weights='distance', p=power)
# Fit the KNeighborsRegressor
knn_model = knn.fit(dats_monthly_jan_2019[['lon', 'lat']], dats_monthly_jan_2019['RR'])
# Predict the value
itrp['RR_2019'] = knn_model.predict(itrp[['lon', 'lat']])
# Read the geojson file
p_jawa = gpd.read_file('data/jawa.geojson')
# Convert predicted data frame into xarray
itrp_xr = xr.DataArray(itrp['RR_2019'].values.reshape(len(grid_lat), len(grid_lon)),
coords={'lat': grid_lat, 'lon': grid_lon},
dims=['lat', 'lon'])
itrp_xr = itrp_xr.rio.write_crs("EPSG:4326")
itrp_xr = itrp_xr.rio.set_spatial_dims("lon", "lat", inplace=True)
# Masking the data with the polygon of Pulau Jawa
itrp_xr_clipped = itrp_xr.rio.clip(p_jawa.geometry, drop=False)
# Plot with matplotlib and cartopy
# Color map customization: red, orange, yellow, green, blue
colors = ['#FF0000', '#FFA500', '#FFFF00', '#00FF00', '#0000FF']
cmap = LinearSegmentedColormap.from_list('mycmap', colors)
# Create a figure and axis
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()}, figsize=(9, 4))
ax.set_extent([106, 114, -8, -6], crs=ccrs.PlateCarree())
# Plot the data
cs = ax.contourf(grid_lon, grid_lat, itrp_xr_clipped.values.reshape(len(grid_lat), len(grid_lon)),
levels=range(0, 700, 25), transform=ccrs.PlateCarree(), cmap=cmap)
fig.colorbar(cs, ax=ax, shrink=0.6, label='Rainfall (mm)')
# Add gridlines just for left and bottom side
ax.set_xticks(np.arange(106, 114, 1), crs=ccrs.PlateCarree())
ax.set_xticklabels([f"{i}°E" for i in np.arange(106, 114, 1)])
ax.set_yticks(np.arange(-8, -5, 1), crs=ccrs.PlateCarree())
ax.set_yticklabels([f"{i}°S" for i in np.arange(8, 5, -1)])
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.add_geometries(p_jawa.geometry, crs=ccrs.PlateCarree(), facecolor='none', edgecolor='black')
ax.title.set_text('Rainfall IDW January 2019')
plt.savefig('idw_py.png', dpi=300, bbox_inches='tight')
plt.show()
Visualisasi Data Iklim (1)
