Plot Raster depth vs value

Specific Domains Geo
1

I’m trying to plot a vertical profile, a 1D plot, of a variable, a multidimensional Raster on a depth-time grid:

using Rasters
using Plots

infile = "sample.nc"
varname = "myvar"

myvar = Raster(infile; name=myvar) # 2D variable
plot(myvar[:, 1]) # vertical profile of timestep 1

This is almost what I want, except that I need the depth axis to be on the vertical axis of the plot, whereas the plot actually produced has the depth on the horizontal axis.

Is there a way to “transpose” the 1D plot?

So far, I’ve found that you can get the desired plot by explicitly supplying the depth axis:

fakedepth = 1:size(myvar,1)
plot(myvar[:,1], fakedepth; yflip=true)

So, I could extract the vertical axis from the netCDF file and then use it to plot the vertical profile . . . but then how does one get hold of the axes in a Raster?

2

It would be best to have a copy of your netcdf to know exactly what is happening here - try to include something we can download.

Normally Rasters.jl will try to plot a Z axis vertically - but I don’t know the names of the dimensions in your data.

To do this manually without Rasters.jl recipes, instead of fakedepth you need lookup(myvar, 1)

To understand lookup a little more and why I don’t say axes: axes is a confusing term because in Julia, axes are a tuple of the indices of each dmension of an array usually e.g. (1:100, 1:200) but for an OffsetArrray might be (0:99, 0:199)

What we want here is not the axes, but their keys/lookups/labels/values whatever is in the dimension variables of your netcdf (what to call this has been bike shedded to death with no winner). In Rasters.jl the dimension variable in the NetCDF is wrapped with a LookupArray for each dimension, as with the regular transforms from GDAL. So I we call them a lookup, for better or worse.

You can get these arrays from a raster with lookup(raster, 1), or better lookup(raster, X), lookup(raster, (X, Y)) etc. They are also what is wrapped in the dimensions returned by dims(raster).

1 Like
3

Thank you for your response! Here is the sample data:

and this is the code:

using Rasters
using Plots
inf = "tmp.nc"
varname = "MYVAR"
@assert(isfile(inf))
myvar = Raster(inf; name=varname)
plot(myvar[:,1])

It plots the depth axis in the horizontal direction.

Normally Rasters.jl will try to plot a Z axis vertically - but I don’t know the names of the dimensions in your data.

It’s named “DEPUV”. It has the attribute axis = "Z", which normally signals that it’s the Z axis. (Is this CF convention?)

I’m not in the position to say what Raster should do, but it would be nice if there are ways to tell Raster which dimension corresponds to Z, etc.

Perhaps it would be too much to ask it to try to “guess” from various clues, like units and dimension names . . . .

In the meanwhile, thank you for the information about lookup! Using it, I’m now able to produce figures I need.

4

Oh, and I’ll be careful about the wording of “axes”.

So, it seems that you struggle with the terminology because the Raster model is more general than netCDF. In netCDF, an “axis” is just a 1D array of real numbers. It’s a “lookup” or “map” from an integer index to the corresponding real number. So, I guess the “lookup” you talk about is much more general than this.

5

On another netCDF file, I tried to do this

v = view(myvar, X(Between(100, 140)), Y(Between(-20,30)))

following the example in Rasters . But, I got

Warning: (X, Y) dims were not found in object

Presumably the reason is the same as for my original problem. In this case, I need a 2D contour plot, but (presumably because of the Warning) plot(v) plots a set of 1D graphs of value vs longitude for all the latitude gridpoints.

This time, I’ve checked the CF convention (NetCDF Climate and Forecast (CF) Metadata Conventions ). It seems to me to be saying that units="degrees_east" is the primary signal that this is latitude. The documentation also says that

the latitude type may be indicated additionally by providing the standard_name attribute with the value latitude, and/or the axis attribute with the value Y.

My latitude has a correct units attribute and axis attribute.

Again, I need a functionality to tell Raster which dimension is X, which is Y, etc.

6

Your dataset doesn’t have X/Y/Z/Time variables or any know conversions to X/Y/Z/Ti.

They are DEPUV and YEARLY1.:

julia> myvar = Raster(inf; name=varname)
105×41 Raster{Union{Missing, Float64},2} MYVAR with dimensions:
  Dim{:DEPUV} Sampled{Float64} Float64[2.5, 7.5, …, 7050.0, 7350.0] ForwardOrdered Explicit Intervals,
  Dim{:YEARLY1} Sampled{DateTime} DateTime[1976-07-02T00:00:00, …, 2016-07-02T00:00:00] ForwardOrdered Explicit Intervals
with missingval: missing
                1976-07-02T00:00:00        1977-07-02T12:00:00  …        2015-07-02T12:00:00        2016-07-02T00:00:00
    2.5   -67460.1                               -78469.4                                  -50027.3                               -34528.7
    7.5   -67093.5                               -77321.2                                  -49179.8                               -34331.5
   12.5   -65447.0                               -74756.2                                  -47710.6                               -32555.5

It does show this right in the REPL :wink:

You can just use the names as keywods in view or getindex/[]

Also IntervalSets.jl is supported now, and .. is more terse than Between:

julia> view(myvar, DEPUV=100..140, YEARLY1=DateTime(1976)..DateTime(1990))
7×14 Raster{Union{Missing, Float64},2} MYVAR with dimensions:
  Dim{:DEPUV} Sampled{Float64} Float64[102.5057795, 107.5404245, …, 128.82229999999998, 134.6500475] ForwardOrdered Explicit Intervals,
  Dim{:YEARLY1} Sampled{DateTime} DateTime[1976-07-02T00:00:00, …, 1989-07-02T12:00:00] ForwardOrdered Explicit Intervals
with missingval: missing
                1976-07-02T00:00:00        1977-07-02T12:00:00  …        1988-07-02T00:00:00        1989-07-02T12:00:00
 102.506  -37600.1                               -36850.3                                  -34487.4                               -40982.9
 107.54   -36761.9                               -35619.4                                  -33620.3                               -40405.4
 112.644  -35316.8                               -33781.6                                  -32292.0                               -39717.7
 117.863  -34392.3                               -32700.4                                  -31898.4                               -39591.6
 123.241  -32614.1                               -30894.6                               …  -31187.4                               -38469.7
 128.822  -30479.4                               -28800.1                                  -29699.8                               -37037.9
 134.65   -28381.0                               -26555.7                                  -27803.2                               -34247.8

You can also write this like:

view(myvar, Dim{:DEPUV}(100..140), Dim{:YEARLY1}(DateTime(1976)..DateTime(1990)))

For plotting, you can put years on the x axis by making yearly a time dimension:

myvar = set(myvar, :YEARLY1 => Ti)
plot(myvar)

To be clear, Rasters is not at all limited to using X/Y/Z/Ti dimension names. They just make a lot of things easier when they are detectable. But you can index with any names. Just read what it says the dimension is in the REPL, or change it with set.

7

It might help to add a keyword for Raster that lets you pass in name changes for the dimensions so you can skip that step. Also detecting units=degrees_east could be done (please write a PR if you have time), although there are more kinds of units for x dimensions than degrees, so its a specific case.

1 Like
8

[Several hours after posting, I fixed some errors in my writing below.] Thank you for your in-depth explanation and tutorial. Here I quickly post this information: The CF convention indicates that the attributes axis="X", axis="Y", etc. are the main means to distinguish the axes. For longitude and latitude, units="degrees_east" still overrides (if it conflicts with the axis attributes). In contrast, non-latitudes and non-longitudes need to rely on the axis attributes.

So, for generality, the axis attribute is the priority to implement before “degrees_east” and “degrees_north”.

Because identification of a coordinate type by its units is complicated by requiring the use of an external software package [UDUNITS] , we provide two optional methods that yield a direct identification. The attribute axis may be attached to a coordinate variable and given one of the values X, Y, Z or T which stand for a longitude, latitude, vertical, or time axis respectively. Alternatively the standard_name attribute may be used for direct identification. But note that these optional attributes are in addition to the required COARDS metadata.

1 Like
9

Ok we can detect axis, it is in your file. The problem with a lot of the CF spec is these are optional attributes and “may be attached” so you cant rely on them.