Open Maschette opened 2 weeks ago
Maschette
commented
2 weeks ago In my for I have added a data frame for what the depth layers are https://github.com/Maschette/bluelink/blob/701237c251e8ac7876035979891a28b0193b0684/data-raw/depth_layers.R
I would like to make either a dataframe that has for each cell what the deepest layer with data is, or just a raster with that. I jsut done know how to do it.
mdsumner
commented
2 weeks ago excellent, nearly all todos I had floating in my head, thanks!!
technically the "windows thing" is that we can only read these successfully on windows via raster package (which reads via ncdf4 package) because the GDAL/terra on CRAN for windows is a bit old. That might change pretty soon so that inconsistency should disappear (we're still in trouble on linux if people have an old GDAL). I didn't really want to workaround all of that, but at least we can read it.
mdsumner
commented
1 week ago ok so depth and time metadata should work now on linux or windows. Oh I changed the function argument to "level", because it's an index not a value.
a <- read_bluelink("2015-06-03", level = 51)
terra::depth(a)
#[1] 4509.18
terra::time(a)
#[1] "2015-06-03 12:00:00 UTCthe WMS is just an image, they normalize it for visualizing as tiles:
terra::rast("WMS:https://thredds.nci.org.au/thredds/wms/gb6/BRAN/BRAN2020/daily/ocean_u_2020_08.nc?service=WMS&version=1.3.0&request=GetCapabilities")
class : SpatRaster
dimensions : 447392431, 1073741824, 3 (nrow, ncol, nlyr)
resolution : 3.352761e-07, 3.352761e-07 (x, y)
extent : -180, 180, -74.95, 75.05 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source : ocean_u_2020_08.nc?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=u&CRS=CRS:84&BBOX=-180.0,-74.95000152638787,180.0,75.05000000050896
colors RGB : 1, 2, 3
names : 84&BBOX=-1~00050896_1, 84&BBOX=-1~00050896_2, 84&BBOX=-1~00050896_3
we can't really use that without other tricks
mdsumner
commented
1 week ago to pull into -180,180 I would do this
temp <- read_bluelink("2023-01-01", varname = "ocean_temp", level = 10)
terra::rotate(temp)
## or more generally
project(temp, rast(ext(-180, 180, -90, 90), res = res(temp), crs = crs(temp))
for subsetting I'll have a think, crop() will work now but won't be as quick on windows
Notes on things to think about/look at.
[ ] This is probably a windows thing but I don't get units in the data on windows but I see it on the main page of the repo.
[ ] Is it also possible to add depth layer as a value there? I don't know enough about SpatRaster to know how to do this. The 'layers' should correspond to the values in
<Dimension name="elevation" units="ncwms:depth" unitSymbol="meters" default="2.5">[ ] Data subsetting: There is this page: https://thredds.nci.org.au/thredds/ncss/grid/gb6/BRAN/BRAN2020/daily/ocean_temp_2014_04.nc/dataset.html which lets you subset the netcdfs down to a certain region to then download them. Not sure if it is useful but at the bottom of it they give you the NCSS request URL which may be useful for if we wanted to subset on the fly with our requests? So if we wanted to subset location down to lon = -180, 0, and lat: -30,-70 the url would be https://thredds.nci.org.au/thredds/ncss/grid/gb6/BRAN/BRAN2020/daily/ocean_temp_2014_04.nc?north=-30.000&west=-180.000&east=0.000&south=-70.000&horizStride=1&time_start=2014-04-01T12:00:00Z&time_end=2014-04-30T12:00:00Z&&&accept=netcdf3
[ ] Projections: I downloaded one of the netcdfs to have a look and they have longitude as -180,180, but when we pull it longitude is 0-360. When I look at the wms data access (eg https://thredds.nci.org.au/thredds/wms/gb6/BRAN/BRAN2020/daily/ocean_u_2020_08.nc?service=WMS&version=1.3.0&request=GetCapabilities) It has the following listed as layers:
and further down has
<BoundingBox CRS="CRS:84" minx="-180.0" maxx="180.0" miny="-74.95000152638787" maxy="75.05000000050896"/>so I dont know if it is possible to pull it in -180,180?