AustralianAntarcticDivision / angstroms

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

polar ROMS #22

Open mdsumner opened 2 years ago

mdsumner commented 2 years ago

waom gridding

f <- " ... ocean_avg_0538-0610_temp.nc"

library(tidync)
x <- tidync(f)
a <- x %>% hyper_filter(ocean_time = index < 2, s_rho = index > 30) %>% hyper_array()

## close: 
#resolution : 2000, 1999.245  (x, y)
#extent     : -3e+06, 3300000, -2699000, 2599000  (xmin, xmax, ymin, ymax)

ex <- raster::extent(-2999000 -1000,  3299000 +1000, -2698000-1000,  2598000+1000)
plot(setExtent(raster(t(a$temp[,ncol(a$temp):1])), ex))

image

mdsumner commented 2 years ago

oh it's 2000 in the y direction, can't see why that would be atm ?

ex <- raster::extent(-2999000 -1000,  3299000 +1000, -2698000-2000,  2598000+2000)
> setExtent(raster(t(a$temp[,ncol(a$temp):1])), ex)
class      : RasterLayer 
dimensions : 2650, 3150, 8347500  (nrow, ncol, ncell)
resolution : 2000, 2000  (x, y)
extent     : -3e+06, 3300000, -2700000, 2600000  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -3.437805, 16.38486  (min, max)
mdsumner commented 2 years ago

VRT for temp (similar will work for u,v,w with some offset/size hacking

(this doesn't work for other ROMS that aren't a regular projected grid)

  tx <- '  <VRTDataset rasterXSize="3150" rasterYSize="2650">
   <GeoTransform>-3e+06, 2000, 0, 2600000, 0, -2000</GeoTransform>
    <SRS>EPSG:3031</SRS>
  <VRTRasterBand dataType="Float32" band="1" blockXSize="450" blockYSize="379">

    <NoDataValue>9.999999933815813e+36</NoDataValue>
    <UnitType>Celsius</UnitType>
    <SimpleSource>
      <SourceFilename relativeToVRT="0">NETCDF:waom2/ocean_avg_0008.nc:temp</SourceFilename>
      <SourceBand>31</SourceBand>
      <SourceProperties RasterXSize="3150" RasterYSize="2650" DataType="Float32" BlockXSize="450" BlockYSize="379" />
      <SrcRect xOff="0" yOff="0" xSize="3150" ySize="2650" />
      <DstRect xOff="0" yOff="0" xSize="3150" ySize="2650" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>'
mdsumner commented 2 years ago

we need some kind of plotting for tidync for this stuff ... but using the warper might actually be better (better than hacky subsetting etc)


get_slice <- function(x, slice = 1L, max = 512) {
  dm <- dim(x)
  idx <- seq_len(prod(dm[1:2]))
  idx <- idx + prod(dm[1:2]) * (slice-1L)
  out <- matrix(x[idx], dm[1:2])
  ix <- seq(1, nrow(out), length.out = max)
  iy <- seq(1, ncol(out), length.out = max)
  out[ix, iy]
}
plot.tidync_data <- function(x, ..., slice = 1L, var = 1L, useRaster = TRUE) {
  image(get_slice(x[[var]]), ..., useRaster = useRaster)
}