weather

#Packages

library(sf)

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(rnaturalearth)

## Warning: package 'rnaturalearth' was built under R version 4.3.3

library(rnaturalearthdata)

## Warning: package 'rnaturalearthdata' was built under R version 4.3.3

## 
## Attaching package: 'rnaturalearthdata'

## The following object is masked from 'package:rnaturalearth':
## 
##     countries110

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(elevatr)

## Warning: package 'elevatr' was built under R version 4.3.3

## elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
## of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
## deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
## will be dropped in future versions, so please plan accordingly.

library(dplyr)
library(chirps)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(nasapower)

## Warning: package 'nasapower' was built under R version 4.3.3

library(terra)

## Warning: package 'terra' was built under R version 4.3.3

## terra 1.8.29

library(stringr)
library(zoo)

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'

## The following object is masked from 'package:terra':
## 
##     time<-

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

#20 random points per each of the 10 domains

library(readr)
library(dplyr)

points_df <- read_csv("InitialDomains/sampled_points.csv")

## Rows: 200 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): point_id
## dbl (3): x, y, domain
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

points_df <- points_df%>%
  rename(lon = x, lat = y, id = point_id)

# Convert to sf object
points_sf <- st_as_sf(points_df, coords = c("lon", "lat"), crs = 4326)
# Get elevation (default source is AWS Terrain Tiles ~90m resolution)
elevation_data <- get_elev_point(points_sf, src = "aws")

## Mosaicing & Projecting

## Note: Elevation units are in meters

# Add altitude to the original data frame
points_with_elev <- points_df %>%
  mutate(elevation = elevation_data$elevation)

points_with_elev <- points_with_elev %>%
  mutate(mz_adaptation = ifelse(elevation >= 1900, "highland", "mid-altitude"))

points_with_elev <- points_with_elev %>%
  mutate(GDD_flo_mz = ifelse(elevation >= 1900, 1053, 1015))

table(points_with_elev$mz_adaptation)

## 
##     highland mid-altitude 
##           64          136

head(points_with_elev)

## # A tibble: 6 × 7
##     lon   lat domain id     elevation mz_adaptation GDD_flo_mz
##   <dbl> <dbl>  <dbl> <chr>      <dbl> <chr>              <dbl>
## 1  29.5 -2.15      1 point1      1823 mid-altitude        1015
## 2  29.7 -2.70      1 point2      1664 mid-altitude        1015
## 3  29.7 -2.39      1 point3      1662 mid-altitude        1015
## 4  29.8 -1.98      1 point4      1597 mid-altitude        1015
## 5  29.9 -2.00      1 point5      1619 mid-altitude        1015
## 6  29.9 -1.96      1 point6      1718 mid-altitude        1015

#Prepare the weather databases with APSIM format, for the last 20 years, from NASA power

dir.create("weather", showWarnings = FALSE)

start_date <- "2005-01-01"
end_date <- "2025-09-01"

for (i in seq_len(nrow(points_with_elev))) {
  point <- points_with_elev[i, ]
  
  try({
    lat <- point$lat
    lon <- point$lon
    id <- point$id

    # === 1. NASA POWER data retrieval ===
    nasapower_data <- get_power(
      community = "AG",
      lonlat = c(lon, lat),
      pars = c("T2M_MAX", "T2M_MIN", "ALLSKY_SFC_SW_DWN", "WS2M", "PRECTOTCORR"),
      dates = c(start_date, end_date),
      temporal_api = "DAILY"
    ) %>%
      rename(
        year = YEAR,
        day = DOY,
        maxt = T2M_MAX,
        mint = T2M_MIN,
        radn = ALLSKY_SFC_SW_DWN,
        wind = WS2M,
        rain = PRECTOTCORR
      ) %>%
      select(year, day, radn, maxt, mint, rain, wind) %>%
      mutate(code = 999999)
    
# === ⚠️ Check for Missing Values ===

      if (anyNA(nasapower_data)) {
      warning(paste("⚠️ Missing values found for point:", id))
      print(nasapower_data %>% filter(if_any(everything(), is.na)) %>% head())
    }
    
    # Fill missing values
nasapower_data <- nasapower_data %>%
  arrange(year, day) %>%
  mutate(
    radn = ifelse(is.na(radn), zoo::rollapplyr(radn, width = 5, FUN = mean, fill = NA, partial = TRUE, align = "right"), radn),
    maxt = ifelse(is.na(maxt), zoo::rollapplyr(maxt, width = 5, FUN = mean, fill = NA, partial = TRUE, align = "right"), maxt),
    mint = ifelse(is.na(mint), zoo::rollapplyr(mint, width = 5, FUN = mean, fill = NA, partial = TRUE, align = "right"), mint),
    wind = ifelse(is.na(wind), zoo::rollapplyr(wind, width = 5, FUN = mean, fill = NA, partial = TRUE, align = "right"), wind),
    rain = ifelse(is.na(rain), 0, rain)
  )
    
    
    # === 2. TAV and AMP ===
    tav <- mean((nasapower_data$maxt + nasapower_data$mint) / 2, na.rm = TRUE)
    amp <- (max((nasapower_data$maxt + nasapower_data$mint)/2, na.rm = TRUE)) -
           (min((nasapower_data$maxt + nasapower_data$mint)/2, na.rm = TRUE))

    # === 3. APSIM format ===
    formatted <- nasapower_data %>%
      transmute(
        year = sprintf("%4d", year),
        day  = sprintf("%3d", day),
        radn = formatC(radn, format = "f", width = 5, digits = 1),
        maxt = formatC(maxt, format = "f", width = 6, digits = 1),
        mint = formatC(mint, format = "f", width = 5, digits = 1),
        rain = formatC(rain, format = "f", width = 5, digits = 1),
        wind = formatC(wind, format = "f", width = 5, digits = 1),
        code = sprintf("%7d", code)
      )

    header <- c(
      "[weather.met.weather]",
      sprintf("latitude = %.4f  (DECIMAL DEGREES)", lat),
      sprintf("longitude = %.4f  (DECIMAL DEGREES)", lon),
      sprintf("tav = %.2f (oC)", tav),
      sprintf("amp = %.2f (oC)", amp),
      "year  day radn  maxt   mint  rain  wind   code
()   () (MJ/m^2) (oC)  (oC)  (mm)  (m/s)    ()"
    )

    file_path <- file.path("weather", paste0(id, ".met"))
    writeLines(header, con = file_path)
    write.table(formatted, file = file_path, append = TRUE, row.names = FALSE,
                col.names = FALSE, quote = FALSE)
    
    message("✅ Saved: ", file_path)
  }, silent = TRUE)
}

## ✅ Saved: weather/point1.met

## ✅ Saved: weather/point2.met

## ✅ Saved: weather/point3.met

## ✅ Saved: weather/point4.met

## ✅ Saved: weather/point5.met

## ✅ Saved: weather/point6.met

## ✅ Saved: weather/point7.met

## ✅ Saved: weather/point8.met

## ✅ Saved: weather/point9.met

## ✅ Saved: weather/point10.met

## ✅ Saved: weather/point11.met

## ✅ Saved: weather/point12.met

## ✅ Saved: weather/point13.met

## ✅ Saved: weather/point14.met

## ✅ Saved: weather/point15.met

## ✅ Saved: weather/point16.met

## ✅ Saved: weather/point17.met

## ✅ Saved: weather/point18.met

## ✅ Saved: weather/point19.met

## ✅ Saved: weather/point20.met

## ✅ Saved: weather/point21.met

## ✅ Saved: weather/point22.met

## ✅ Saved: weather/point23.met

## ✅ Saved: weather/point24.met

## ✅ Saved: weather/point25.met

## ✅ Saved: weather/point26.met

## ✅ Saved: weather/point27.met

## ✅ Saved: weather/point28.met

## ✅ Saved: weather/point29.met

## ✅ Saved: weather/point30.met

## ✅ Saved: weather/point31.met

## ✅ Saved: weather/point32.met

## ✅ Saved: weather/point33.met

## ✅ Saved: weather/point34.met

## ✅ Saved: weather/point35.met

## ✅ Saved: weather/point36.met

## ✅ Saved: weather/point37.met

## ✅ Saved: weather/point38.met

## ✅ Saved: weather/point39.met

## ✅ Saved: weather/point40.met

## ✅ Saved: weather/point41.met

## ✅ Saved: weather/point42.met

## ✅ Saved: weather/point43.met

## ✅ Saved: weather/point44.met

## ✅ Saved: weather/point45.met

## ✅ Saved: weather/point46.met

## ✅ Saved: weather/point47.met

## ✅ Saved: weather/point48.met

## ✅ Saved: weather/point49.met

## ✅ Saved: weather/point50.met

## ✅ Saved: weather/point51.met

## ✅ Saved: weather/point52.met

## ✅ Saved: weather/point53.met

## ✅ Saved: weather/point54.met

## ✅ Saved: weather/point55.met

## ✅ Saved: weather/point56.met

## ✅ Saved: weather/point57.met

## ✅ Saved: weather/point58.met

## ✅ Saved: weather/point59.met

## ✅ Saved: weather/point60.met

## ✅ Saved: weather/point61.met

## ✅ Saved: weather/point62.met

## ✅ Saved: weather/point63.met

## ✅ Saved: weather/point64.met

## ✅ Saved: weather/point65.met

## ✅ Saved: weather/point66.met

## ✅ Saved: weather/point67.met

## ✅ Saved: weather/point68.met

## ✅ Saved: weather/point69.met

## ✅ Saved: weather/point70.met

## ✅ Saved: weather/point71.met

## ✅ Saved: weather/point72.met

## ✅ Saved: weather/point73.met

## ✅ Saved: weather/point74.met

## ✅ Saved: weather/point75.met

## ✅ Saved: weather/point76.met

## ✅ Saved: weather/point77.met

## ✅ Saved: weather/point78.met

## ✅ Saved: weather/point79.met

## ✅ Saved: weather/point80.met

## ✅ Saved: weather/point81.met

## ✅ Saved: weather/point82.met

## ✅ Saved: weather/point83.met

## ✅ Saved: weather/point84.met

## ✅ Saved: weather/point85.met

## ✅ Saved: weather/point86.met

## ✅ Saved: weather/point87.met

## ✅ Saved: weather/point88.met

## ✅ Saved: weather/point89.met

## ✅ Saved: weather/point90.met

## ✅ Saved: weather/point91.met

## ✅ Saved: weather/point92.met

## ✅ Saved: weather/point93.met

## ✅ Saved: weather/point94.met

## ✅ Saved: weather/point95.met

## ✅ Saved: weather/point96.met

## ✅ Saved: weather/point97.met

## ✅ Saved: weather/point98.met

## ✅ Saved: weather/point99.met

## ✅ Saved: weather/point100.met

## ✅ Saved: weather/point101.met

## ✅ Saved: weather/point102.met

## ✅ Saved: weather/point103.met

## ✅ Saved: weather/point104.met

## ✅ Saved: weather/point105.met

## ✅ Saved: weather/point106.met

## ✅ Saved: weather/point107.met

## ✅ Saved: weather/point108.met

## ✅ Saved: weather/point109.met

## ✅ Saved: weather/point110.met

## ✅ Saved: weather/point111.met

## ✅ Saved: weather/point112.met

## ✅ Saved: weather/point113.met

## ✅ Saved: weather/point114.met

## ✅ Saved: weather/point115.met

## ✅ Saved: weather/point116.met

## ✅ Saved: weather/point117.met

## ✅ Saved: weather/point118.met

## ✅ Saved: weather/point119.met

## ✅ Saved: weather/point120.met

## ✅ Saved: weather/point121.met

## ✅ Saved: weather/point122.met

## ✅ Saved: weather/point123.met

## ✅ Saved: weather/point124.met

## ✅ Saved: weather/point125.met

## ✅ Saved: weather/point126.met

## ✅ Saved: weather/point127.met

## ✅ Saved: weather/point128.met

## ✅ Saved: weather/point129.met

## ✅ Saved: weather/point130.met

## ✅ Saved: weather/point131.met

## ✅ Saved: weather/point132.met

## ✅ Saved: weather/point133.met

## ✅ Saved: weather/point134.met

## ✅ Saved: weather/point135.met

## ✅ Saved: weather/point136.met

## ✅ Saved: weather/point137.met

## ✅ Saved: weather/point138.met

## ✅ Saved: weather/point139.met

## ✅ Saved: weather/point140.met

## ✅ Saved: weather/point141.met

## ✅ Saved: weather/point142.met

## ✅ Saved: weather/point143.met

## ✅ Saved: weather/point144.met

## ✅ Saved: weather/point145.met

## ✅ Saved: weather/point146.met

## ✅ Saved: weather/point147.met

## ✅ Saved: weather/point148.met

## ✅ Saved: weather/point149.met

## ✅ Saved: weather/point150.met

## ✅ Saved: weather/point151.met

## ✅ Saved: weather/point152.met

## ✅ Saved: weather/point153.met

## ✅ Saved: weather/point154.met

## ✅ Saved: weather/point155.met

## ✅ Saved: weather/point156.met

## ✅ Saved: weather/point157.met

## ✅ Saved: weather/point158.met

## ✅ Saved: weather/point159.met

## ✅ Saved: weather/point160.met

## ✅ Saved: weather/point161.met

## ✅ Saved: weather/point162.met

## ✅ Saved: weather/point163.met

## ✅ Saved: weather/point164.met

## ✅ Saved: weather/point165.met

## ✅ Saved: weather/point166.met

## ✅ Saved: weather/point167.met

## ✅ Saved: weather/point168.met

## ✅ Saved: weather/point169.met

## ✅ Saved: weather/point170.met

## ✅ Saved: weather/point171.met

## ✅ Saved: weather/point172.met

## ✅ Saved: weather/point173.met

## ✅ Saved: weather/point174.met

## ✅ Saved: weather/point175.met

## ✅ Saved: weather/point176.met

## ✅ Saved: weather/point177.met

## ✅ Saved: weather/point178.met

## ✅ Saved: weather/point179.met

## ✅ Saved: weather/point180.met

## ✅ Saved: weather/point181.met

## ✅ Saved: weather/point182.met

## ✅ Saved: weather/point183.met

## ✅ Saved: weather/point184.met

## ✅ Saved: weather/point185.met

## ✅ Saved: weather/point186.met

## ✅ Saved: weather/point187.met

## ✅ Saved: weather/point188.met

## ✅ Saved: weather/point189.met

## ✅ Saved: weather/point190.met

## ✅ Saved: weather/point191.met

## ✅ Saved: weather/point192.met

## ✅ Saved: weather/point193.met

## ✅ Saved: weather/point194.met

## ✅ Saved: weather/point195.met

## ✅ Saved: weather/point196.met

## ✅ Saved: weather/point197.met

## ✅ Saved: weather/point198.met

## ✅ Saved: weather/point199.met

## ✅ Saved: weather/point200.met