Error with ERA5 data in skewT plot for altitudes > 2000 meters

abigmo commented 2 years ago

I am trying to produce skew T plot from ERA5 data. I am mainly interested in 2021 and 2022 data, since in rawinsonde are missing.

Using some ERA5 data (sample from here )

I have this error with the latest thundeR.

era5 <- readRDS("sample_era5.RDS") 

r.earth <- 6.3781E+6

era5 %>%
    dplyr::filter(time == time[20] & level >= 750) %>% ## select one time
    plyr::arrange(desc(level)) %>% ## make sure it's arranged for increasing height
    mutate("wd" = mod(180 + 180/pi*atan2(v, u), 360)) %>% ## ERA5 has wind "blowing to" ## https://confluence.ecmwf.int/pages/viewpage.action?pageId=133262398
    mutate("wd" = mod(wd + 180, 360)) %>% ## convert wind direction to blowing from, since thundeR wants "blowing from"
    mutate("ws" = ( sqrt(u^2 + v^2) ) / 0.51444) %>% ## wind speed in knots (1 knots = 0.51444 m/s)
    mutate("dewpoint" = dewpoint.calc(p = level, temp = t, rh = r) - 273.15) %>% ## compute dew point with custom formula
    mutate("height" = z * r.earth / (9.80665 * r.earth - z)) %>% ## from geopot to height with correction for decrease in g with altitude
    with(., sounding_save(filename = "test.png", pressure = level, altitude = height, temp = t - 273.15, dpt = dewpoint, wd = wd, ws = ws, accuracy = 2))

Error in if (output$pressure[output$altitude - output$altitude[1] == 2000] >  :
  argument is of length zero

Basically it throws an error any time the altitude is larger than 2000 meters: is is a problem from my script or of thundeR?

Two extra feedbacks for the developers:

the sounding data needs to be sorted for increasing altitude in order to be processed correctly by the package: it's OK, but it's better to say it explicitly in the manual
it looks like that the wind direction in input needs to be indicated as "blowing from": is that correct?

Any help is appreciated. And also the option to select 2021 and 2022 in rawinsonde would be great.

Thanks, AB

PS: version.string R version 3.6.3 (2020-02-29) platform x86_64-pc-linux-gnu
arch x86_64

bczernecki commented 2 years ago

Hi @abigmo ,

I have briefly looked at the sample dataset and it seems to be working for some of the selected test cases. However, I wasn't able to compute wind direction and dew point using your code as it seems to be inherited from some external packages or other chunks of your code. Could you send a link to a complete gist that will be fully reproducible so we can investigate it with more details?

Thanks, Bartek

bczernecki commented 2 years ago

Two extra feedbacks for the developers:

the sounding data needs to be sorted for increasing altitude in order to be processed correctly by the package: it's OK,

it looks like that the wind direction in input needs to be indicated as "blowing from": is that correct?

Yes, both informations are correct. Entirely agree with adding these information to user guides for sake of clarity. Thanks for this hint!

abigmo commented 2 years ago

@bczernecki thanks for the prompt feedback :-) I have fiddled some more with the units of rh and temp, just to make sure I did not made any extra mistake. I actually have found some, but there might be more...but I still have this mistake with the altitude > 2000 m.

I am copy-pasting everything again below, including the extra code for the dew.point.

Any help is appreciated. Thanks!

Alessandro

## -------------------------------------------------------------------------------------------------------------------------------------------------------
##  File Title:   fun_dew_point.R
##  Description:  Calculate dew point from rh an T_Air
##  Autor:        Christian Brida  
##                Institute for Alpine Environment
##  Data:         26/07/2017
##  Version:      1.0
## ------------------------------------------------------------------------------------------------------------------------------------------------------

##  Method 4: GEOtop formula https://github.com/geotopmodel/geotop/blob/master/src/geotop/meteo.cc

## Jan 2021: minor modifications by Alessandro Bigi - University of Modena and Reggio Emilia

## Saturated vapour pressure 
## 
## @param temp Air Temperature (°C)
## @param p atm pressure (hPa)
##
dewpoint.calc <- function(p = NULL, temp.celsius = NULL, rh.percent = NULL){

    SatVapPressure <- function(p = p, temp = temp.celsius){
        a <- 6.1121 * (1.0007 + p * 3.46E-6)
        b <- 17.502
        c <- 240.97
        e <- a * exp((b * temp) / (c + temp))
        return(e)
    }

    ## Temperature from saturated water vapour 
    ## 
    ## @param SVP Saturated vapour pressure
    ## @param P Air pressure  
    ##
    TfromSatVapPressure <- function(SVP, p = p){
        a <- 6.1121 * (1.0007 + p * 3.46e-6)
        b <- 17.502
        c <- 240.97
        temp.dew <- c * log(SVP / a) / (b - log(SVP / a))
        return (temp.dew)
    }

    ## compute saturated vapor pressure
    e <- SatVapPressure(temp = temp.celsius, p = p)

    ## rh[rh > 100] <- 100
    ## rh[rh < 0] <- 0
    ## compute the temp corresponding to saturated vapor pressure
    svp <- e * rh.percent / 100

    dew.celsius <- TfromSatVapPressure(SVP = svp, p = p)

    return(dew.celsius)
}

era5 <- readRDS("sample_era5.RDS") 

r.earth <- 6.3781E+6
library("thunder")

era5 %>%
    dplyr::filter(time == time[20] & level >= 750) %>% ## select one time
    plyr::arrange(desc(level)) %>% ## make sure it's arranged for increasing height
    mutate("wd" = magrittr::mod(180 + 180/pi*atan2(v, u), 360)) %>% ## wind blowing to! ## https://confluence.ecmwf.int/pages/viewpage.action?pageId=133262398
    mutate("wd" = magrittr::mod(wd + 180, 360)) %>% ## convert wind direction to blowing from, since thunder wants "blowing from"
    mutate("ws" = ( sqrt(u^2 + v^2) ) / 0.51444) %>% ## wind speed in knots (1 knots = 0.51444 m/s)
    mutate("dewpoint.celsius" = dewpoint.calc(p = level, temp.celsius = t - 273.15, rh.percent = r)) %>% ## compute dew point in celsius with custom formula
    mutate("height" = z * r.earth / (9.80665 * r.earth - z)) %>% ## from geopot to height with correction for decrease in g with altitude
    with(., sounding_save(filename = "test.png", pressure = level, altitude = height, temp = t - 273.15, dpt = dewpoint.celsius, wd = wd, ws = ws, accuracy = 2))

bczernecki commented 2 years ago

I hope I have succesfully patched all problems. The most recent changes were committed into the devel branch. The reproducible example is given below. Hope that works!

Thanks again for the feedback!

# install devel version:
devtools::install_github("bczernecki/thundeR@devel")
library(thunder)
library(tidyverse)

download.file("https://box.ing.unimo.it/owncloud/index.php/s/mJeL6zwUEGFx81o/download",
              mode = "wb",
              "~/era.rds")

era5 <- readRDS("~/era.rds") 
r.earth <- 6.3781E+6
dewpoint.calc <- function(p = NULL, temp.celsius = NULL, rh.percent = NULL){

  SatVapPressure <- function(p = p, temp = temp.celsius){
    a <- 6.1121 * (1.0007 + p * 3.46E-6)
    b <- 17.502
    c <- 240.97
    e <- a * exp((b * temp) / (c + temp))
    return(e)
  }

  ## Temperature from saturated water vapour 
  ## 
  ## @param SVP Saturated vapour pressure
  ## @param P Air pressure  
  ##
  TfromSatVapPressure <- function(SVP, p = p){
    a <- 6.1121 * (1.0007 + p * 3.46e-6)
    b <- 17.502
    c <- 240.97
    temp.dew <- c * log(SVP / a) / (b - log(SVP / a))
    return (temp.dew)
  }

  ## compute saturated vapor pressure
  e <- SatVapPressure(temp = temp.celsius, p = p)

  ## rh[rh > 100] <- 100
  ## rh[rh < 0] <- 0
  ## compute the temp corresponding to saturated vapor pressure
  svp <- e * rh.percent / 100

  dew.celsius <- TfromSatVapPressure(SVP = svp, p = p)

  return(dew.celsius)
}

df = era5 %>%
  #dplyr::filter(time == time[20] & level >= 750) %>% ## select one time
  dplyr::filter(time == time[20] & level >= 100) %>% ## select everything up to 100hPa
  plyr::arrange(desc(level)) %>% ## make sure it's arranged for increasing height
  mutate("wd" = magrittr::mod(180 + 180/pi*atan2(v, u), 360)) %>% ## wind blowing to! ## https://confluence.ecmwf.int/pages/viewpage.action?pageId=133262398
  mutate("wd" = magrittr::mod(wd + 180, 360)) %>% ## convert wind direction to blowing from, since thunder wants "blowing from"
  mutate("ws" = ( sqrt(u^2 + v^2) ) / 0.51444) %>% ## wind speed in knots (1 knots = 0.51444 m/s)
  mutate("dewpoint.celsius" = dewpoint.calc(p = level, temp.celsius = t - 273.15, rh.percent = r)) %>% ## compute dew point in celsius with custom formula
  mutate("height" = z * r.earth / (9.80665 * r.earth - z)) ## from geopot to height with correction for decrease in g with altitude

with(df, sounding_plot(pressure = level, altitude = height, temp = t-273.15, dpt = dewpoint.celsius, wd = wd, ws = ws, accuracy = 2))

# optionally you can clip to any other height
df2 = dplyr::filter(df, level<= 800)

with(df2, sounding_plot(pressure = level, altitude = height, temp = t-273.15, dpt = dewpoint.celsius, wd = wd, ws = ws, accuracy = 2))

abigmo commented 2 years ago

Thanks for the prompt fix. It works like a charm :+1:

bczernecki / thundeR

Error with ERA5 data in skewT plot for altitudes > 2000 meters #38