bài tập

Make something

hahjagđhfd

library(ncdf4) library(raster) library(sp) library(dplyr) library(sf) library(leaflet) library(openxlsx) # Set working directory setwd("c:/users/daniel/downloads") # Open the NetCDF files d18O <- nc_open("LGMR_d18Op_climo.nc") SAT <- nc_open("LGMR_SAT_climo.nc") # Extract dimensions and variables d18O_obsdatadates <- d18O$dim$age$vals lon <- ncvar_get(d18O, "lon") lat <- ncvar_get(d18O, "lat") d18O_time <- ncvar_get(d18O, "age") d18O_tmp_array <- ncvar_get(d18O, "d18Op") SAT_tmp_array <- ncvar_get(SAT, "sat") SAT_fillvalue <- ncatt_get(SAT, "sat", "_FillValue")$value # Fill value for missing data d18O_fillvalue <- ncatt_get(d18O, "d18Op", "_FillValue")$value # Create a data frame of coordinates from the NetCDF grid lonlat <- as.matrix(expand.grid(lon, lat)) # Define SP data frame SP <- data.frame( lon = as.numeric(c("47.55", "45.366", "41.848", "42.076")), lat = as.numeric(c("19.45", "22.873", "20.675", "20.896")) ) # Convert SP to SpatialPoints coordinates(SP) <- ~lat + lon proj4string(SP) <- CRS("+init=epsg:4326") # Prepare a data frame to store results results <- data.frame() # Loop through each time slice for (t in 1:length(d18O_obsdatadates)) { # Extract the slice for the current time step for both datasets d18O_nc_slice <- d18O_tmp_array[, , t] d18O_nc_slice[d18O_nc_slice == d18O_fillvalue] <- NA # Replace fill values with NA SAT_nc_slice <- SAT_tmp_array[, , t] SAT_nc_slice[SAT_nc_slice == SAT_fillvalue] <- NA # Replace fill values with NA # Convert the d18O slice to a data frame d18O_nc_vec <- as.vector(d18O_nc_slice) d18O_VP <- data.frame(cbind(lonlat, d18O_nc_vec)) names(d18O_VP) <- c("lon", "lat", "d18Op") d18O_VP <- na.omit(d18O_VP) # Convert the SAT slice to a data frame SAT_nc_vec <- as.vector(SAT_nc_slice) SAT_VP <- data.frame(cbind(lonlat, SAT_nc_vec)) names(SAT_VP) <- c("lon", "lat", "sat") SAT_VP <- na.omit(SAT_VP) # Combine d18O and SAT data frames combined_VP <- merge(d18O_VP, SAT_VP, by = c("lon", "lat")) # Convert combined_VP to SpatialPointsDataFrame coordinates(combined_VP) <- ~lon + lat proj4string(combined_VP) <- CRS("+init=epsg:4326") # Find the nearest neighbor for each point in SP nearest_indices <- apply(coordinates(SP), 1, function(pt) { distances <- spDistsN1(as.matrix(coordinates(combined_VP)), pt, longlat = TRUE) which.min(distances) }) # Extract the nearest values nearest_d18Op <- combined_VP@data[nearest_indices, "d18Op"] nearest_sat <- combined_VP@data[nearest_indices, "sat"] # Combine results with time and SP coordinates temp_results <- data.frame( lon = coordinates(SP)[, 1], lat = coordinates(SP)[, 2], time = d18O_obsdatadates[t], d18Op = nearest_d18Op, sat = nearest_sat ) # Append to the main results data frame results <- rbind(results, temp_results) } # Close the NetCDF files nc_close(d18O) nc_close(SAT) # Save results to a CSV file write.csv(results, "nearest_values.csv", row.names = FALSE) # Create a map plot for the first age slice d18O_first_age_slice <- d18O_tmp_array[, , 1] d18O_first_age_slice[d18O_first_age_slice == d18O_fillvalue] <- NA d18O_first_age_df <- data.frame(cbind(lonlat, as.vector(d18O_first_age_slice))) colnames(d18O_first_age_df) <- c("lon", "lat", "d18Op") d18O_first_age_df <- na.omit(d18O_first_age_df) # Use sf package for easier handling of spatial data # Convert d18O_first_age_df to sf object d18O_first_age_sf <- st_as_sf(d18O_first_age_df, coords = c("lon", "lat"), crs = 4326) sat_first_age_slice <- SAT_tmp_array[, , 1] sat_first_age_slice[sat_first_age_slice == SAT_fillvalue] <- NA sat_first_age_df <- data.frame(cbind(lonlat, as.vector(sat_first_age_slice))) colnames(sat_first_age_df) <- c("lon", "lat", "d18Op") sat_first_age_df <- na.omit(sat_first_age_df) # Use sf package for easier handling of spatial data # Convert d18O_first_age_df to sf object sat_first_age_sf <- st_as_sf(sat_first_age_df, coords = c("lon", "lat"), crs = 4326) # Convert SP to a data frame SP_df <- data.frame( lon = SP@coords[, 1], lat = SP@coords[, 2], extracted_d18Op = results$d18Op[results$time == d18O_obsdatadates[1]], # Add extracted d18Op values for the first age extracted_sat = results$sat[results$time == d18O_obsdatadates[1]] # Add extracted SAT values for the first age ) # Convert SP_df to sf object SP_sf <- st_as_sf(SP_df, coords = c("lon", "lat"), crs = 4326) # Create a base plot using leaflet leaflet() %>% addTiles() %>% addCircleMarkers( data = d18O_first_age_sf %>% st_coordinates() %>% as.data.frame(), lng = ~X, lat = ~Y, radius = 4, color = "blue", fillOpacity = 0.6, label = ~paste("d18Op:", round(d18O_first_age_df$d18Op, 2), round(sat_first_age_df$d18Op, 2)) ) %>% addCircleMarkers( data = SP_sf %>% st_coordinates() %>% as.data.frame(), lng = ~X, lat = ~Y, radius = 6, color = "red", fillOpacity = 1, label = ~paste("SP Point d18Op:", round(SP_df$extracted_d18Op, 2), "SAT:", round(SP_df$extracted_sat, 2)) ) # Save results to a CSV file write.xlsx(results, "netCDF_extract_Zoli.xlsx", rowNames = FALSE)

Cities by day and night