R Notebook

# Increase download timeout
options(timeout = 300)

# Load libraries
library(terra)

## terra 1.9.50

library(geodata)

# Create folders
if (!dir.exists("Data")) dir.create("Data")
if (!dir.exists("Output")) dir.create("Output")

# Download Bangladesh National Boundary (Level 0)
bd_boundary <- gadm(country = "BGD", level = 0, path = "Data")

# Plot boundary
plot(bd_boundary, main = "Bangladesh National Boundary")

# Download 30 arc-second elevation data
bd_elev <- elevation_30s(country = "BGD", path = "Data")

# Crop elevation data using boundary
bd_elev_crop <- crop(bd_elev, bd_boundary)

# Mask elevation data to country shape
bd_elev_mask <- mask(bd_elev_crop, bd_boundary)

# Plot elevation map
plot(bd_elev_mask, main = "Elevation Map of Bangladesh (m)")

# Check resolution
res(bd_elev_mask)

## [1] 0.008333333 0.008333333

# Check Coordinate Reference System (CRS)
crs(bd_elev_mask)

## [1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        MEMBER[\"World Geodetic System 1984 (G2296)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

# Save processed raster to local drive as GeoTIFF
writeRaster(bd_elev_mask, filename = "Output/elevation_bd.tif", overwrite = TRUE)

# 1. Load required library
library(terra)

# 2. Load the saved GeoTIFF raster from local drive
elev_bd <- rast("Output/elevation_bd.tif")

# 3. Calculate minimum, maximum, and mean elevation values using global()
elev_min  <- global(elev_bd, fun = "min", na.rm = TRUE)
elev_max  <- global(elev_bd, fun = "max", na.rm = TRUE)
elev_mean <- global(elev_bd, fun = "mean", na.rm = TRUE)

# Print the calculated values
print("--- MINIMUM ELEVATION ---")

## [1] "--- MINIMUM ELEVATION ---"

# 1. Increase elevation values of the raster by 10 meters
elev_plus_10 <- elev_bd + 10

# 2. Example: Convert a raw temperature raster by dividing values by 10
# Assuming 'raw_temp' is a raw temperature raster (e.g., WorldClim temp * 10)
# temp_converted <- raw_temp / 10

# Demonstrating map algebra division on our raster dataset:
elev_divided <- elev_bd / 10

# Display the modified elevation raster (+10 meters)
plot(elev_plus_10, main = "Elevation Increased by 10 Meters")

# 1. Define the reclassification matrix:
# From -> To -> Become (Class)
# 0 to 10    -> Class 1
# 10 to 50   -> Class 2
# 50 to Inf  -> Class 3
reclass_matrix <- matrix(c(0, 10, 1,
                           10, 50, 2,
                           50, Inf, 3), 
                         ncol = 3, 
                         byrow = TRUE)

# 2. Apply reclassification using classify()
elev_classified <- classify(elev_bd, reclass_matrix)

# 3. Plot the classified raster map
plot(elev_classified, main = "Reclassified Elevation Classes (1, 2, 3)")

# 1. Load ggplot2 library
library(ggplot2)

# 2. Load the built-in mtcars dataset
data(mtcars)

# 3. Create base ggplot object and add a histogram layer for mpg
a <- ggplot(mtcars, aes(x = mpg))
a + geom_histogram(binwidth = 5, fill = "skyblue", color = "black") +
  labs(title = "Histogram of Miles Per Gallon (mpg)", x = "Miles Per Gallon (mpg)", y = "Count")

# 1. Scatter plot of mpg (X-axis) vs disp (Y-axis) with linear trend line
ggplot(mtcars, aes(x = mpg, y = disp)) +
  geom_point(color = "blue", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  labs(title = "Scatter Plot of mpg vs disp with Linear Trend Line",
       x = "Miles Per Gallon (mpg)",
       y = "Engine Displacement (disp)")

## `geom_smooth()` using formula = 'y ~ x'

# 1. Load plotly library for interactivity
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# 2. Create extended scatter plot with aesthetics mapping
p <- ggplot(mtcars, aes(x = mpg, y = disp, color = factor(carb), shape = factor(carb), size = wt)) +
  geom_point() +
  labs(title = "Interactive Scatter Plot of mpg vs disp",
       x = "Miles Per Gallon (mpg)",
       y = "Engine Displacement (disp)",
       color = "Carburetors",
       shape = "Carburetors",
       size = "Weight (wt)")

# 3. Convert ggplot object to interactive plot using ggplotly()
ggplotly(p)