Hello there! I am having troubles in this script to plot SST maps and help to make this one run again is very much apreciated... I will upload the.nc and script to anyone thats can solve this issue
link to .nc file
https://mega.nz/file/bAIXWZSK#oRbuiKCYqvKrGxYWDUckUwrAVT9bVej7zMssLsxGYsE
# Leitura do diretório onde os arquivos estão
setwd("please change to your directory")
# Livrarias necessárias
require(ncdf4)
require(reshape2)
library(ggplot2)
library(sf)
library(leaflet)
library(dplyr)
library(leaflet.extras)
library(rayshader)
require(tidyverse)
require(ncdump)
library(oce)
library(lubridate)
library(metR)
require(metR)
require(tidyverse)
require(lubridate)
require(oce)
require(ocedata)
require(sf)
library(ggspatial)
metadata = ncdump::NetCDF("SST_17_04_2022.nc")
# Ler o arquivo NCDF
nc <- nc_open("SST_17_04_2022.nc")
nc
## Extrair componente espacial
lat <- ncvar_get(nc, "lat")
lon <- ncvar_get(nc, "lon")
## Extrair componente temporal
time = ncvar_get(nc, "time")
# convert time original (to) to julian
to = insol::JDymd(year = 1950, month = 1, day = 1)
# add the original time to the extracted time
jd = to+time
#convert the julian day to gregorian calender
date = insol::JD(jd, inverse = TRUE)
## sea surface temperature
sst<- ncvar_get(nc, "analysed_sst")
## atribuindo dimnames conforme lon, lat e tempo
dimnames(sst)[[1]] <- as.character(lon)
dimnames(sst)[[2]] <- as.character(lat)
dimnames(sst)[[3]] <- as.character(date)
str(sst)
## melt para reestruturação dos dados
sst1 <- melt(sst)
head(sst1,50)
summary(sst1)
colnames(sst1)<-c("lon","lat","sst")
summary(sst1)
# Selecionando a área
db2<-subset(sst1,sst1$lat<=5&lat>=0&lon<=-39&lon>=-46)
summary(db2)
## convert drifter observation into simple features
drifter.split<-db2
drifter.split.sf = drifter.split %>%
st_as_sf(coords = c("lon", "lat"))
## divide the Atlantic ocean region into grids
drifter.grid = drifter.split.sf %>%
st_make_grid(n = c(70,60))%>%
st_sf()
#drifter.split.sf.se = drifter.split.sf%>% filter(season=="SE")
drifter.gridded = drifter.grid %>%
mutate(id = 1:n(), contained = lapply(st_contains(st_sf(geometry),drifter.split.sf),identity),
obs = sapply(contained, length),
sst = sapply(contained, function(x) {median(drifter.split.sf[x,]$sst, na.rm = TRUE)}))
drifter.gridded = drifter.gridded %>% select(obs,sst) %>% na.omit()
## obtain the centroid coordinates from the grid as table
coordinates = drifter.gridded %>%
st_centroid() %>%
st_coordinates() %>%
as_tibble() %>%
rename(x = X, y = Y)
## remove the geometry from the simple feature of gridded drifter dataset
st_geometry(drifter.gridded) = NULL
## stitch together the extracted coordinates and drifter information int a single table for SE monsoon season
current.gridded.se = coordinates %>%
bind_cols(drifter.gridded) %>%
mutate(season = "SE")
## bind the gridded table for SE and NE
## Note that similar NE follow similar procedure, hence not shown in the post
drifter.current.gridded = current.gridded.se %>%
bind_rows(current.gridded.se)
## select grids for SE season only
drf.se = drifter.current.gridded %>%
filter(season == "SE")
## interpolate the U component
sst.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$sst)
## obtain dimension that determine the width (ncol) and length (nrow) for tranforming wide into long format table
dimension = data.frame(lon = sst.se$xg, sst.se$zg) %>% dim()
## make a sst data table from interpolated matrix
sst.tb = data.frame(lon = sst.se$xg,
sst.se$zg) %>%
gather(key = "lata", value = "sst", 2:dimension[2]) %>%
mutate(lat = rep(sst.se$yg, each = dimension[1])) %>%
select(lon,lat, sst) %>% as_tibble()
## Error is below this line(Also Needs to convert sst from Kelvin to Celsius)
## stitch now the V component intot the U data table and compute the velocity
uv.se = sst.tb %>%
bind_cols(sst.tb %>% select(sst))
# Criando breaks para legendas
breaks = seq(-1, 32, by = 0.5)
# Planejamento
# Operação
#longitude<-c(-31.5386,-30.3426,-30.4137,-31.0345)
#latitude<-c( 1.5748,1.8921,2.5548,2.7328)
#lo<-as.data.frame(longitude)
#la<-as.data.frame(latitude)
#plan<-cbind.data.frame(longitude,latitude)
##########################################################################################################
# Mapas com apenas contornos de SST e a SST
p<-ggplot() +
# Criar o mapa de SST
metR::geom_contour_fill(data = uv.se, aes(x = lon, y = lat, z = sst),
na.fill = TRUE, bins = 250) +
# Criar o mapa de contorno (isolinhas) de SST
metR::geom_contour2(data = uv.se, aes(x = lon, y = lat, z = sst),color="black",size=0.5)+
# overlay dos contornos (isolinhas) de SST no mapa
metR::geom_text_contour(data = uv.se, aes(x = lon, y = lat, z = sst),
parse = TRUE, check_overlap = TRUE, size = 3)+
# Adicionar o mapa
geom_sf(data = spData::world, fill = "black", col = 1)+
# limitar o mapa as coordenadas
coord_sf(xlim = c(-39,-46), ylim = c(5,0))+
# Adicionar pontos de interesse (no caso pesca)
#geom_point(data=plan,aes(x =longitude, y =latitude),color="blue", size=5)+
# Adicionar a escala de valores das variaveis( limits = c(-0.2,0.2),
# breaks =seq(-0.2,0.2,.05)) e cores (colours = oce::oce.colorsJet(120))
scale_fill_gradientn(name = "SST\n(ºC)",colours = c("#caacc3","#0000c4","#0039f9","#0355f7","#00facd","#03fe8d","#01ff2a",
"#41ff01","#bdfe00","#f3b704","#dd6431","#d12d30","#9f0e30","#452109",
"#271102"),breaks=breaks,na.value = NA)+
theme_bw()+
# configuração de letras e fontes do mapa
theme(legend.position = "right",
legend.key.height = unit(1.4, "cm"),
legend.background = element_blank(),
axis.text = element_text(size = 12, colour = 1))+
# Adicionar titulo, subtitulo e caption ao mapa
labs(
title = "GOES 16 12:00 N5 15/05/2020",
subtitle = "Sea Surface Temperature (ºC)",
caption = "Creation: SATFISHING SARVS",
x = expression(paste("Longitude ", "(", degree, ")")),
y = expression(paste("Latitude ", "(", degree, ")"))) +
#scale_x_continuous(breaks = seq(-31.5, -25.5, .5)) +
#scale_y_continuous(breaks = seq(0.5, 4.5, .5)) +
#Adicionar grades de acordo com a latitude e longitude
theme(
panel.background = element_rect(fill = NA),
panel.grid.major = element_line(colour = "grey80"), # Adiciona cores as linhas ex. black...
panel.ontop = TRUE
)
# Salvar em BMP
bmp(filename = "SST_17_04_2022_SSTCONT.bmp",
res=250,width = 3000, height = 1800)
print(p)
dev.off()
#############################################################################################################
#Plotar apenas SST
p<-ggplot() +
geom_raster(data = uv.se, aes(x = lon, y = lat, fill = sst), interpolate = TRUE)+
geom_sf(data = spData::world, fill = "black", col = 1)+
labs(title = "FISHING GROUNDS N5 15/05/2020",
subtitle = "Sea Surface Temperature (ºC)",
caption = "Creation: SATFISHING SARVS") +
coord_sf(xlim = c(-39,-46), ylim = c(5,0))+
#geom_point(data=plan,aes(x =longitude, y =latitude),color="black", size=12)+
scale_fill_gradientn(name = "SST\n(ºC)",colours =c("#caacc3","#0000c4","#0039f9","#0355f7","#00facd","#03fe8d","#01ff2a",
"#41ff01","#bdfe00","#f3b704","#dd6431","#d12d30","#9f0e30","#452109",
"#271102"), breaks=breaks,na.value = NA)+
theme_bw()+
theme(legend.position = "right",
legend.key.height = unit(1.4, "cm"),
legend.background = element_blank(),
axis.text = element_text(size = 12, colour = 1))+
labs(x = "", y = "")
# Salvar em BMP
bmp(filename = "SST_17_04_2022.bmp",
res=250,width = 3000, height = 1800)
print(p)
dev.off()
Erro: Aesthetics must be either length 1 or the same as the data (2091): z
Run `rlang::last_error()` to see where the error occurred.
> dev.off()
null device
1