harphub / harpIO

IO functions for HARP
https://harphub.github.io/harpIO/
Other
6 stars 16 forks source link

Calculate wind speed/direction from u,v components #37

Closed meteorolog90 closed 11 months ago

meteorolog90 commented 3 years ago

Sorry if I bother you with such questions, I don't have much experience with the R language.

Does harp have a function to calculate the wind speed and direction from the u and v wind components ?I tried to use functions from meteogrid/R/geowind.R, but I didn't know how to work with it.

What i did: -read u,v components from SQlite:

u10 <- read_point_forecast (

  start_date = 2020082000,
  end_date   = 2020082100,
  fcst_model = "ALAEF",
  fcst_type  = "eps", 
  parameter  = "u10m",
  file_path  = "/work/users/ext005/sql/forecast",

)

v10 <- read_point_forecast (

  start_date = 2020082500,
  end_date   = 2020082600,
  fcst_model = "ALAEF",
  fcst_type  = "eps",
  parameter  = "v10m",
  file_path  = "/work/users/ext005/sql/forecast",

)

and then I used the functions

> wind.dirspeed(u10,v10)
Error in u^2 : non-numeric argument to binary operator
> geowind(u10,v10)
Error in geowind(u10, v10) : u must be a geofield!

I guess the problem is in object u10, v10, but I don't know how to change them.

andrew-MET commented 3 years ago

There are a couple of issues here.

The meteogrid functions are to be applied on gridded data as they require projection information for wind direction, so need to be applied before interpolation in the read_forecast() step. Unfortunately this isn't built into the grib reading yet - hopefully @adeckmyn can address that soon.

harp actually doesn't know what do with wind direction when verifying at the moment as it won't take account of the cyclic nature of wind direction. I have opened an issue (https://github.com/harphub/harpPoint/issues/22) for that, but it might not be addressed until next year due to lack of harp budget for this year.

However, you can get wind speed from the data you have - it will require an update of harpPoint to be able to do this.

library(harp)
library(dplyr) # for rename, mutate, select & inner_join
library(purrr) # for map2

u10 <- read_point_forecast (

  start_date = 2020082000,
  end_date   = 2020082100,
  fcst_model = "ALAEF",
  fcst_type  = "eps", 
  parameter  = "u10m",
  file_path  = "/work/users/ext005/sql/forecast",

)

v10 <- read_point_forecast (

  start_date = 2020082500,
  end_date   = 2020082600,
  fcst_model = "ALAEF",
  fcst_type  = "eps",
  parameter  = "v10m",
  file_path  = "/work/users/ext005/sql/forecast",

)

# Get all the members in a single column for easy manipulation
# and remove the parameter column so it's not used in the join

u10 <- gather_members(u10) %>% 
  rename(u = forecast) %>% 
  select(-parameter)

v10 <- gather_members(v10) %>% 
  rename(v = forecast) %>% 
  select(-parameter)

# get wind speed by joining and calculating the square root of the sum of the squares
# making sure that the output is of class "harp_fcst" and return to each member having its own column
# You can optionally add a new parameter column, but it's not actually used for anything.

wind_speed <- map2(u10, v10, inner_join) %>% 
  structure(class = "harp_fcst") %>% 
  mutate(forecast = sqrt(u ^ 2 + v ^ 2), parameter = "S10m") %>% 
  select(-u, -v) %>% 
  spread_members()

# Maybe get rid of u10 and v10 to save memory
rm(u10, v10)

You can then do verification of wind speed.

adeckmyn commented 3 years ago

The read_fa() function in harpIO already has a way to get wind speed or direction from (u,v) components. I may be able to implement the same system in read_grib. Note by the way that calculating wind speed as sqrt(u^2+v^2) doesn't take account of the map factor, so it is only approximately correct.

adeckmyn commented 3 years ago

With grib files, an extra difficulty is that the code must then also ascertain that the U and V components are along the model grid or whether they were already rotated to lat/lon axis. I believe there is a key for this, but accessing it needs a direct call to Rgrib2 which will make the code very complex. Just having it as a function argument (like I did for FA: fa_rotate=TRUE), is easier. For speed it doesn't matter so much, but for direction obviously it does.

andrew-MET commented 3 years ago

I think the likelihood of this being fully addressed this year is quite small, so in the meantime I would suggest post processing the grib files before harp to get wind speed and direction using gl:

https://hirlam.org/trac/wiki/HarmonieSystemDocumentation/40h1.1/PostPP/gl

If you don't have access see:

http://hirlam.org/index.php/register

meteorolog90 commented 3 years ago

@andrew-MET thanks a lot. I was thinking in this direction, but what i said, i have not a lot experience in R programming, so i didint know how to deal with this problem. What i did:

join two objects: wind <- u10$ALAEF %>% left_join(v10$ALAEF,by=c("SID","fcdate","fcst_cycle","leadtime","units","validdate","z"))

drop, change, replace parameter name:

wind <- wind %>% select(-parameter.y)
wind <- wind %>% rename(parameter=parameter.x)
wind <- mutate(wind,parameter=ifelse(parameter=="u10m","S10m"))

calculate wind speed: wind <- wind %>% mutate(ALAEF_mbr000=(ALAEF_mbr000.x)^2+(ALAEF_mbr000.y)^2)

drop columns:

wind <- wind %>% select(-ALAEF_mbr000.y)
wind <- wind %>% select(-ALAEF_mbr000.x)

with this approach, structure for object wind was change, and lost the name for tibble -> got error when i wanna to save verification in to rds format. (mname not exist).

> str(wind)
tibble [4,624 x 12] (S3: tbl_df/tbl/data.frame)
 $ SID         : int [1:4624] 11060 11080 11800 11803 11805 11806 11810 11812 11813 11815 ...
 $ fcdate      : int [1:4624] 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 ...
 $ fcst_cycle  : chr [1:4624] "00" "00" "00" "00" ...
 $ leadtime    : int [1:4624] 0 0 0 0 0 0 0 0 0 0 ...
 $ parameter   : chr [1:4624] "S10m" "S10m" "S10m" "S10m" ...
 $ units       : chr [1:4624] "m/s" "m/s" "m/s" "m/s" ...
 $ validdate   : int [1:4624] 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 1598313600 ...
 $ z           : int [1:4624] 10 10 10 10 10 10 10 10 10 10 ...
 $ ALAEF_mbr000: num [1:4624] 2.62 0.674 1.45 0.862 0.46 ...
 $ ALAEF_mbr001: num [1:4624] 1.808 1.882 0.14 0.177 0.656 ...
 $ ALAEF_mbr002: num [1:4624] 2.8049 0.9464 0.9141 0.0646 0.4635 ...
 $ ALAEF_mbr003: num [1:4624] 0.976 0.954 4.551 0.296 0.555 ...

That is why I gave up this approach and looked for other solutions (meteogrid / R / geowind.R). When i was fail with these geowind functions either, I gave up and wrote you a message.

With you code, object wind_speed have right structure(that's what I was looking for):

> wind_speed
* ALAEF
# A tibble: 4,624 x 12
     SID fcdate fcst_cycle leadtime units validdate     z parameter ALAEF_mbr000
   <int>  <int> <chr>         <int> <chr>     <int> <int> <chr>            <dbl>
 1 11060 1.60e9 00                0 m/s      1.60e9    10 S10m             1.62
 2 11060 1.60e9 00                3 m/s      1.60e9    10 S10m             1.72
 3 11060 1.60e9 00                6 m/s      1.60e9    10 S10m             1.32
 4 11060 1.60e9 00                9 m/s      1.60e9    10 S10m             2.18
 5 11060 1.60e9 00               12 m/s      1.60e9    10 S10m             3.37
 6 11060 1.60e9 00               15 m/s      1.60e9    10 S10m             2.45
 7 11060 1.60e9 00               18 m/s      1.60e9    10 S10m             0.897
 8 11060 1.60e9 00               21 m/s      1.60e9    10 S10m             1.98
 9 11060 1.60e9 00               24 m/s      1.60e9    10 S10m             2.00
10 11060 1.60e9 00               27 m/s      1.60e9    10 S10m             2.01
# ... with 4,614 more rows, and 3 more variables: ALAEF_mbr001 <dbl>,
#   ALAEF_mbr002 <dbl>, ALAEF_mbr003 <dbl>

Now saving point verification scores working!