# Plot Raster depth vs value

**URL:** <https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150>\
**Category:** Geo\
**Created:** [October 23, 2022, 3:58pm UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150 "2022-10-23T15:58:39Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [October 23, 2022, 3:58pm UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/1 "2022-10-23T15:58:39Z")

</div>

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

```julia
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:

```julia
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`?

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 23, 2022, 4:56pm UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/2 "2022-10-23T16:56:32Z")

</div>

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)`.

---

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [October 24, 2022, 5:40am UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/3 "2022-10-24T05:40:41Z")

</div>

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

> **[tmp.nc](https://drive.google.com/file/d/1uwRIqcbNNvU-Be_jVCf3ZU-QxniXaHnU/view?usp=sharing)**
>
> Google Drive file.

and this is the code:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [October 24, 2022, 5:58am UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/4 "2022-10-24T05:58:46Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [October 29, 2022, 2:33pm UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/5 "2022-10-29T14:33:55Z")

</div>

On another netCDF file, I tried to do this

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

```

following the example in [Rasters](https://rafaqz.github.io/Rasters.jl/dev/) . But, I got

```julia
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](https://cfconventions.org/cf-conventions/cf-conventions.html#latitude-axis-ex) ). 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.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 29, 2022, 2:43pm UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/6 "2022-10-29T14:43:45Z")

</div>

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
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 😉

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
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:

```julia
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:

```julia
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`.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 29, 2022, 3:26pm UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/7 "2022-10-29T15:26:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [November 5, 2022, 5:14am UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/8 "2022-11-05T05:14:49Z")

</div>

[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.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [November 5, 2022, 10:44am UTC](https://discourse.julialang.org/t/plot-raster-depth-vs-value/89150/9 "2022-11-05T10:44:27Z")

</div>

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.