AustralianAntarcticDivision / angstroms

Tools for model grid data (with raster in R)
http://australianantarcticdivision.github.io/angstroms/
4 stars 1 forks source link

How are the different extents dealt with for u and v? #10

Open RMORSEcode opened 6 years ago

RMORSEcode commented 6 years ago

Hi Mike,

I am getting NA values in some fields for either u or v (but not both) and I think it is because the extents of u and v vectors are different (they are offset by 1 pixel in either u or v direction for the respective variable), so when a cell is chosen based on the face_z_index from the bgm file, the corresponding u or v may not actually be in that cell. I am using 24-hour averaged roms output, so there should not be any native NA values for just one of the vectors.

ru <- set_indextent(brick(roms_file, varname = "u", lvar = 4, level = level, ncdf=T)) rv <- set_indextent(brick(roms_file, varname = "v", lvar = 4, level = level, ncdf=T))

extent(ru) class : Extent xmin : 0 xmax : 241 ymin : 0 ymax : 106

extent(rv) class : Extent xmin : 0 xmax : 242 ymin : 0 ymax : 105

face_z_index$v <- extract_at_level(rv, rename(face_z_index, level = roms_level, cell = cell))

head(face_z_index)

atlantis_level roms_level face cell u v

1 3 1 face0 10162 NA -0.025098816 2 3 2 face0 10162 NA -0.016623283 3 2 3 face0 10162 NA 0.002793798 4 2 4 face0 10162 NA 0.029781619 5 2 5 face0 10162 NA 0.029100351 6 2 6 face0 10162 NA 0.026689434
RMORSEcode commented 6 years ago

image

I guess the simple answer is to make them the same - but what does this actually do to the data???

ncol(ru) [1] 241 ncol(rv) [1] 242 ncol(rv)=241 nrow(ru) [1] 106 nrow(rv) [1] 105 nrow(ru)=105 dim(ru) [1] 105 241 40 dim(rv) [1] 105 241 40

This should allow for comparisons... but I get an error about the data not being in memory:

face_z_index$u <- extract_at_level(ru, rename(face_z_index, level = roms_level, cell = cell)) Error in .readCells(x, cells, 1) : no data on disk or in memory

mdsumner commented 6 years ago

You can't do that, it's throwing the data away.

I think this means we have to build the index independently for u and v, which means we need to map the faces into the different coordinate spaces as well

## UNTESTED
roms_ll_u <- romscoords(roms_file, transpose = TRUE, spatial = c("lon_u", "lat_u"))
roms_ll_v <- romscoords(roms_file, transpose = TRUE, spatial = c("lon_v", "lat_v"))
roms_face_u <- romsmap(project_to(faceSpatial(bgm), "+init=epsg:4326"), roms_ll_u)
roms_face_v <- romsmap(project_to(faceSpatial(bgm), "+init=epsg:4326"), roms_ll_v)

ind_face_u <- cellFromLine(romsdata(roms_file, "u"), roms_face_u)
ind_face_v <- cellFromLine(romsdata(roms_file, "v"), roms_face_v)
face_roms_index <- tibble(face = roms_face$label[rep(seq_len(nrow(roms_face)), lengths(ind_face))], 
                          cell_u = unlist(ind_face_u),                           cell_v = unlist(ind_face_v))

and so on all the way through. I'd just set this aside when I was working on it before, but it's clearer now how to go about it. I don't think I can spend any time on this in the next two weeks though, I'm afraid. I'll see what I can do thoug

RMORSEcode commented 6 years ago

Great, I will work with that and let you know what happens, thanks for replying so quickly. The salinity and temperature are calculated at the internal rho points, which have lat_rho and lon_rho, so those would be used for anything that uses box methods, rather than the face methods.

Another option to normalize the grid would be to interpolate between upper and lower v values for a given cell and left and right u values, respectively on a layer by layer basis. We would have to deal with issues where one of the values may not exist (so either drop cell or use singular value).