Upgrade to Pro — share decks privately, control downloads, hide ads and more …

harp_spatial_inca.pdf

Christoph
October 17, 2019

 harp_spatial_inca.pdf

harp spatial inca observations

Christoph

October 17, 2019
Tweet

Other Decks in Education

Transcript

  1. How to get your spatial observation to harp Many institutes

    have their own tool for spatial observation - analysis, and it is most probably not coming in a format that harp understands right a way (hdf5, grib, fa, fatar). So what to do?
  2. How to get your spatial observation to harp Many institutes

    have their own tool for spatial observation - analysis, and it is most probably not coming in a format that harp understands right a way (hdf5, grib, fa, fatar). So what to do? meteogrid is a good starting point to help you constructing this spatial observation.
  3. library( readr ) library( tidyverse ) # file is gzipped

    inca <- read_delim( incafile, delim = " ", col_names = FALSE ) %>% dplyr::mutate_all( as.numeric ) How do we read the data
  4. Inca output is actually starting in the bottom left, going

    along the top left line by line - BUT It is written to the ascii file in a really weird way: • Starting bottom left, • writing line by line • but only 10 cols (instead of 701) inca <- c( t( inca ))[!is.na(c( t( inca )))] %>% matrix( nrow = 701 ) Re-arranging needed
  5. R is indexing rows first, second is columns so •

    Transposing the matrix • Make a vector out of it -> get all inca data in one stream • Remove all NA’s • Turn it into a matrix again inca <- c( t( inca ))[!is.na(c( t( inca )))] %>% matrix( nrow = 701 ) Re-arranging needed
  6. Tell R what map - projection the matrix has Proj4string

    from https://epsg.io/ +proj=lcc +lat_0=47.5N +lon_0=13.3333 +lat_1=46N +lat_2=49N +ellps=bessel +x_0=400000 +y_0=400000 701 x 401 domain Projection summary: proj= lcc NE = ( 17.74378 , 49.39726 ) SW = ( 8.444457 , 45.77268 ) Now make a geofield out of it
  7. domainINCA <- structure( list( projection = list( proj = "lcc",

    lat_0 = 47.5, lon_0 = 13.3333, lat_1 = 46, lat_2 = 49, ellps = "bessel", x_0 = 400000, y_0 = 400000), nx = 701, ny = 401, SW = c(8.444457, 45.772682), NE = c(17.74378, 49.397259) ), class = "geodomain") Now make a geofield out of it
  8. To be read by the library( harp ) environment( read_inca

    ) <- environment( read_grid ) Put it into the environment