plot WRF: Cross sections

Based on this NCL

library(eixport)
library(raster)
#> Loading required package: sp
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1
library(cptcity)
library(sf)
library(vein)
library(ggplot2)

Reading Temperature and crop coast lines for our study area

wrfo <- "/media/sergio/ext41/wrfo_sebrazil/wrfout_d02_2014-11-03_000000"

t2 <- wrf_get(wrfo, "T2", as_raster = T)
#> Warning in wrf_get(wrfo, "T2", as_raster = T): The option as_raster will be
#> deprecated, see eixport::wrf_raster()
t2 <- t2 -273.15

Find colour palette for temperature

find_cpt("temperature")
#> [1] "arendal_temperature"    "idv_temperature"        "jjg_misc_temperature"  
#> [4] "kst_03_red_temperature"

Let us create a line between c(-46.5,-23.85) and c(-46.35, -23.95)

m <- cbind(c(-46.5, -46.35),  # xini xend
           c(-23.85, -23.95)) # yini yend
cross = st_linestring(m)
(cross <- st_sfc(cross, crs = 4326))
#> Geometry set for 1 feature 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -46.5 ymin: -23.95 xmax: -46.35 ymax: -23.85
#> Geodetic CRS:  WGS 84
#> LINESTRING (-46.5 -23.85, -46.35 -23.95)
plot(t2$T2_2014.11.03_00.00.00, 
     main = "Temperature using plot", 
     col = cpt("arendal_temperature"))
plot(cross, add = T)

Now, define several lines

m2 <- cbind(c(-46.05, -46.36),  # xini xend
            c(-23.85, -23.95))   # yini yend

Define a helper function

points_extract <- function(m, sta) {
  cross = st_linestring(m)
  cross <- st_sfc(cross, crs = 4326)
  t2s <- st_as_sf(sta)
  lt <- st_intersection(t2s, cross) 
  geo <- st_geometry(lt)
  lt <- st_set_geometry(lt, NULL)
  na <- names(lt)
  lt$id <- 1:nrow(lt)
  dx <- vein::wide_to_long(df = lt, 
                           column_with_data = na, 
                           column_fixed = "id")
  stf <- st_sf(dx, geometry = geo)
  lt <- st_centroid(stf)
  lt <- cbind(lt, st_coordinates(lt))
  
  return(lt)
}

sta = st_as_stars(t2)
names(sta) <- "temperature"
df <- points_extract(m, sta = sta)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
#> Warning in data.frame(..., check.names = FALSE): row names were found from a
#> short variable and have been discarded
#> Warning in st_centroid.sf(stf): st_centroid assumes attributes are constant over
#> geometries of x

Let us check the data

head(df)
#> Simple feature collection with 6 features and 5 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -46.49466 ymin: -23.90841 xmax: -46.41243 ymax: -23.85356
#> CRS:           +proj=longlat +datum=WGS84 +no_defs
#>         V1 V2                     V3         X         Y
#> 1 17.93331  1 T2_2014.11.03_00.00.00 -46.49466 -23.85356
#> 2 18.83734  2 T2_2014.11.03_00.00.00 -46.48186 -23.86211
#> 3 20.90524  3 T2_2014.11.03_00.00.00 -46.46182 -23.87548
#> 4 22.21670  4 T2_2014.11.03_00.00.00 -46.44715 -23.88527
#> 5 22.05804  5 T2_2014.11.03_00.00.00 -46.43040 -23.89643
#> 6 22.49484  6 T2_2014.11.03_00.00.00 -46.41243 -23.90841
#>                      geometry
#> 1 POINT (-46.49466 -23.85356)
#> 2 POINT (-46.48186 -23.86211)
#> 3 POINT (-46.46182 -23.87548)
#> 4 POINT (-46.44715 -23.88527)
#> 5  POINT (-46.4304 -23.89643)
#> 6 POINT (-46.41243 -23.90841)

Add time variable, select and plot

df$time <- as.POSIXct(df$V3, format = "T2_%Y.%m.%d_%H.%M.%S")
df$hour <- strftime(df$time, "%H")
ggplot(df[df$hour %in% "15", ], 
       aes(x = X, y = V1, colour = as.factor(time))) +
  labs(y =expression(paste("Temperature [",degree,"C]")),
       x = expression(paste("Longitude [",degree,"]")))+
  geom_line() +
  theme_bw()+
  theme(legend.title = element_blank())

Sergio Ibarra-Espinosa

2021-07-17