# TrueColor raster satellite image

**URL:** <https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401>\
**Category:** Geo\
**Tags:** gmt, gdal, rasters\
**Created:** [March 9, 2024, 5:22pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401 "2024-03-09T17:22:07Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 9, 2024, 5:22pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/1 "2024-03-09T17:22:07Z")

</div>

I’m trying to compose a TrueColor image from three raster channels a VIIRS granule stored in a netCDF4 file. I’d like to do this in pure Julia but haven’t figured it out yet. I can read the data using NCDatasets.jl, but can’t figure out the composing of the RGB image so I can plot with GeoMakie.jl using Rasters.jl or the like in color constructed from three input channels.

The alternative I have tried is to use GMT.jl and RemoteS.jl. The example on the RemoteS.jl page is elegant, however everything I’ve tried to read in the data using gmtread() has been unsuccessful. My best attempt is below…

So I’m stuck. Does anyone have any suggestions for getting something like this working. I’m open to any and all solutions.

Many thanks

```julia
# following https://github.com/GenericMappingTools/RemoteS.jl
using RemoteS, GMT

Ib = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname=:observation_data/M03, grd=true);
Ig = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname=:observation_data/M04, grd=true);
Ir = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname=:observation_data/M05, grd=true);
Irgb = truecolor(Ir, Ig, Ib);
imshow(Irgb)

```

---

<div class="post-metadata">

**Author:** ![Fliks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fliks/32/2494_2.png) [@Fliks](https://discourse.julialang.org/u/Fliks)\
**Post date:** [March 9, 2024, 5:58pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/2 "2024-03-09T17:58:38Z")

</div>

You can use Rasters and NCDatasets to load the data and then construct an RGB image from that.

```julia
using Rasters
using NCDatasets
using Colors

ras = Raster("pathtodata.nc")
rgb = RGB.(view(ras, Band=1)./255, view(ras, Band=2)./255, view(ras, Band=3)./255)

```

Here you can select the Bands or however your dimension is named with the normal Selectors from DimensionalData.  
If your data contains missing values you can overload the RGB of three missing inputs, but that seems to be a bit dangerous.

```julia
Colors.RGB(::Missing, ::Missing, ::Missing) = RGB(0,0,0)

```

Now you should be able to plot the data using GLMakie or GeoMakie.

```julia
using GLMakie
plot(rgb)

```

I am not sure how good this would work when your data is too large for your RAM.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 9, 2024, 6:20pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/3 "2024-03-09T18:20:25Z")

</div>

> [@haffner](#):
>
> however everything I’ve tried to read in the data using gmtread() has been unsuccessful.

I’m not sure that a `Symbol` like `:observation_data/M03` containing a slash is supposed to work. Note that in the above example the `truecolor` function is expecting a UInt16 image (don’t remember if other types are accepted too). If `gmtread` is failing maybe `gdalread` can save your day. Do you have a pointer to the `"VJ102MOD.A2022017.0218.021.2022017084708.nc"` file so that I can try it?

---

<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:** [March 9, 2024, 8:57pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/4 "2024-03-09T20:57:58Z")

</div>

Should we add a helper method to Rasters.jl to do this? seems like a common use case.

Like `to_rgb(rast; dims=Band(1:3)) = ...`

Then you could select difference bands or even a different dimension. It could handle setting missing values to alpha transparency.

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 9, 2024, 9:32pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/5 "2024-03-09T21:32:14Z")

</div>

@joa-quim Thanks for the reply. Storing netCDF4 datasets in groups is fairly common Here’s a link to a file. The values should be readily converted to UInt16 if necessary.

[https://acd-ext.gsfc.nasa.gov/anonftp/toms/tmp/dhaffner/VJ102MOD.A2022017.0218.021.2022017084708.nc](https://acd-ext.gsfc.nasa.gov/anonftp/toms/tmp/dhaffner/VJ102MOD.A2022017.0218.021.2022017084708.nc)

---

<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:** [March 9, 2024, 9:40pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/6 "2024-03-09T21:40:15Z")

</div>

What joa-quim is refering to is this is not possible syntax:

```julia
:observation_data/M05

```

You may mean

```julia
"observation_data/M05"

```

But did you try @Fliks example above? you may need to tweak it for you dimension names of you netcfd, but otherwise that should work.

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 9, 2024, 9:48pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/7 "2024-03-09T21:48:14Z")

</div>

Having trouble loading the data with this simple Raster call. The bands are stored in separate netCDF datasets inside a parent group, “observation\_data”, like this:

```julia
group: observation_data {
  variables:
  	ushort M01(number_of_lines, number_of_pixels) ;
  	ushort M02(number_of_lines, number_of_pixels) ;
  	ushort M03(number_of_lines, number_of_pixels) ;
  	ushort M04(number_of_lines, number_of_pixels) ;
  	ushort M05(number_of_lines, number_of_pixels) ;
  	ushort M06(number_of_lines, number_of_pixels) ;
	...

```

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 9, 2024, 9:59pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/8 "2024-03-09T21:59:03Z")

</div>

I tried a syntax with the colon prefix because that was documented in the readwrite subprogram of GMT.jl

I tried this instead and was unable to get the read to work.

```julia
Ib = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname="observation_data/M03", grd=true);

```

Reading the data is not an issue though, since I can read using NCDatasets with this

```julia
ds = Dataset("../data/viirs/n20/VJ102MOD.A2022017.0218.021.2022017084708.nc", "r")
M3 = ds.group["observation_data"]["M03"]
M4 = ds.group["observation_data"]["M04"]
M5 = ds.group["observation_data"]["M05"]

```

I appreciate your help - once I’ve learned the basics of working with this particular type of satellite data, I’ll draft an intro tutorial for others.

---

<div class="post-metadata">

**Author:** ![Fliks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fliks/32/2494_2.png) [@Fliks](https://discourse.julialang.org/u/Fliks)\
**Post date:** [March 9, 2024, 10:56pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/9 "2024-03-09T22:56:44Z")

</div>

I tried to open the data with YAXArrays and Rasters and they both don’t understand the structure of this particular dataset.

Is this a common way to structure a NetCDF file as a collection of subgroups?

I opened an issue at Rasters.jl

> <https://github.com/rafaqz/Rasters.jl/issues/605>
>
> We can't open the NetCDF file from this link as as Raster even though you can op…en it with NCDatasets. 
> https://acd-ext.gsfc.nasa.gov/anonftp/toms/tmp/dhaffner/VJ102MOD.A2022017.0218.021.2022017084708.nc
> 
> The Netcdf file is nested and the actual data is organised in groups at a deeper level.
> One can construct a Raster, by opening the data with NCDatasets and using one of the subgroups but that is cumbersome.
> 
> See this discourse discussion;
> https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/8

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 9, 2024, 11:29pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/10 "2024-03-09T23:29:57Z")

</div>

Well, I got a moderate success.

```julia
I3 = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname="observation_data/M03", gdal=true);
Warning 1: dimension #1 (number_of_pixels) is not a Longitude/X dimension.
Warning 1: dimension #0 (number_of_lines) is not a Latitude/Y dimension.

I4 = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname="observation_data/M04", gdal=true);
Warning 1: dimension #1 (number_of_pixels) is not a Longitude/X dimension.
Warning 1: dimension #0 (number_of_lines) is not a Latitude/Y dimension.

I5 = gmtread("VJ102MOD.A2022017.0218.021.2022017084708.nc", varname="observation_data/M05", gdal=true);
Warning 1: dimension #1 (number_of_pixels) is not a Longitude/X dimension.
Warning 1: dimension #0 (number_of_lines) is not a Latitude/Y dimension.

Irgb = truecolor(I3, I4, I5);

```

Note, one can use `gdalread` and it’s even a bit simpler syntax (but used code is the same)

```julia
I3 = gdalread("NETCDF:VJ102MOD.A2022017.0218.021.2022017084708.nc:/observation_data/M03");

```

Now, this is doing an automatic scaling to UInt8 and I need to check if it’s doing automatic histogram stretch or not.

The result is not very exciting but although it seems to be a gray image, it’s in fact a RGB color one, but one where layers have very similar value and hence it looks gray. The white lines are in the data.

```julia
viz(Irgb)

```

 ![GMTjl_j-or8](https://global.discourse-cdn.com/julialang/original/3X/b/9/b9e0a968180767b09e65efa6c881c41420c18738.jpeg)

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 10, 2024, 1:33pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/11 "2024-03-10T13:33:23Z")

</div>

It’s common to organize datasets in groups in the larger netCDF files used in the Earth Observation community (EU Copernicus and NASA). The instrument data and derived science data products are often complex - the groups help to organize things.

Thanks for opening an issue. I can provide additional example files if desired (e.g. from TROPOMI).

---

<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:** [March 10, 2024, 1:40pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/12 "2024-03-10T13:40:18Z")

</div>

[https://cfconventions.org/cf-conventions/cf-conventions.html#groups](https://cfconventions.org/cf-conventions/cf-conventions.html#groups)

The real question is how can we represent these with DimensionalData.jl (so Rasters.jl/YAXArrays.jl etc). We don’t have a concept of nestedness. It would be like a `DimStack` of `DimStack`. Maybe that will work without much hassle.

Or is it better to flatten and combine the names like `:observation_data_MO3`

For now we can also just add a `group` keyword to get a Raster/RasterStack for a specific group, without the nesting

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 10, 2024, 2:55pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/13 "2024-03-10T14:55:30Z")

</div>

Adding the group keyword seems good if you can also specify the exact dataset. Reading all datasets in a group may work in some cases but not all. Many datasets are 2-d so they could stack easily, but some datasets are 3-d (native raster stacks) while others may be 1-d. Users specify which to read that are conformable, or read them each manually and create their own raster stacks as they see fit maybe? Perhaps this is moving away from the way things are done in Rasters.jl though…

PS I don’t think there any need to capture the nestedness of the original datasets in the internal (to Julia) representation of the data.

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 10, 2024, 3:23pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/14 "2024-03-10T15:23:08Z")

</div>

That looks similar to my image using my old methods in IDL, though I get a color image and it’s reversed left-to-right. I’ll try to get color out of your code as a starting point. I appreciate it.

I think the lines are expected because of how the instrument samples data in the flight direction (S-N) to accommodate larger fields of view towards the edges of the scan.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/f/df8af5a8cf55a6faa9f8dce55db1791abc4e6bfd.jpeg)

---

<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:** [March 10, 2024, 3:23pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/15 "2024-03-10T15:23:19Z")

</div>

Rasters/DimensionalData can mix layer dimensions, they dont all have to be 2d.

Another question: if there are groups, does loading the first group with info on how to load the other groups with a keyword sound like a good default?

---

<div class="post-metadata">

**Author:** ![haffner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/haffner/32/35871_2.png) [@haffner](https://discourse.julialang.org/u/haffner)\
**Post date:** [March 10, 2024, 3:34pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/16 "2024-03-10T15:34:16Z")

</div>

That sounds like a good idea for files where the first subgroup contains only datasets, but some groups may have a mix of datasets and groups at a second level. Not sure how to handle that effectively by default. Users who know the layout of the file can pick and choose though.

---

<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:** [March 10, 2024, 4:08pm UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/17 "2024-03-10T16:08:50Z")

</div>

Oh right they actually nest that deeply in practice?

Damn. We probably need a new type for this.

---

<div class="post-metadata">

**Author:** ![Fliks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fliks/32/2494_2.png) [@Fliks](https://discourse.julialang.org/u/Fliks)\
**Post date:** [March 11, 2024, 9:07am UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/18 "2024-03-11T09:07:02Z")

</div>

As a intermediate solution you can construct a Raster from the subgroup that you are interested in.  
And then you can use the scaleminmax function from ImageCore to scale your image and construct a colorview for it.

```julia
using Rasters
using NCDatasets
using Colors, ImageCore

julia> nc = NCDataset("VJ102MOD.A2022017.0218.021.2022017084708.nc")
julia> obs = nc.group["observation_data"]
julia> M3 = Raster(obs, path, :M03, missingval = obs["M03"][1,1])
julia> M4 = Raster(obs, path, :M04, missingval = obs["M04"][1,1])
julia> M5 = Raster(obs, path, :M05, missingval = obs["M05"][1,1])

```

Now we compute the extrema of the raster that we then use for scaling the colorview linearly

```julia
julia> m3stats = extrema(skipmissing(M3))
julia> m4stats = extrema(skipmissing(M4))
julia> m5stats = extrema(skipmissing(M5))

julia> rgbscale = colorview(RGB, scaleminmax(m3stats...).(M3), scaleminmax(m4stats...).(M4), scaleminmax(m5stats...).(M5))

```

I am not sure whether the coloring of the image is as expected. You might want to play around with the minima and maxima setting and rather use the quantile function to get some not as extreme values.

Constructing the Raster objects that way feels a bit hacky and it would be good when the missing value could be constructed directly from the metadata of the Arrays.

---

<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:** [March 11, 2024, 9:15am UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/19 "2024-03-11T09:15:40Z")

</div>

To be clear thats not documented API sytax 😉

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 11, 2024, 10:56am UTC](https://discourse.julialang.org/t/truecolor-raster-satellite-image/111401/20 "2024-03-11T10:56:34Z")

</div>

> [@haffner](#):
>
> That looks similar to my image using my old methods in IDL,

Can you describe what that method does?

The problem with this data is that the bands are UInt16 arrays but they have scale **two** factors, and I’m emphasizing the plural because each band has a `scale_factor` and a `radiance_scale_factor`. Which one of the two should be applied? Next issue is that after applying a scale factor we no longer have a UInt16 but instead a Float array. And how then do the necessary histogram stretch that is necessary to make a color image really look like a color image?

You can do this to see the bands metadata (only showing a part of it)

```julia
gdalinfo "NETCDF:VJ102MOD.A2022017.0218.021.2022017084708.nc:/observation_data/M03"

...
Band 1 Block=3200x16 Type=UInt16, ColorInterp=Undefined
  NoData Value=65535
  Offset: 0, Scale:1.999175765377e-05
  Metadata:
    add_offset=0
    flag_meanings=Missing_EV Bowtie_Deleted Cal_Fail
    flag_values={65532,65533,65534}
    long_name=M-band 03 observations at pixel locations
    NETCDF_VARNAME=M03
    radiance_add_offset=0
    radiance_scale_factor=0.013002133
    radiance_units=Watts/meter^2/steradian/micrometer
    scale_factor=1.9991758e-05
    valid_max=65527
    valid_min=0
    _FillValue=65535

```

If I do not apply the scale factor (which GMT.jl does by default when if finds one in a bands metadata) and skip the `_FillValue`’s and compute the band’s histogram plus the automatic limits for the histogram stretching, I get this (not very different for the M4 & M5 bands).

 ![GMTjl_j](https://global.discourse-cdn.com/julialang/original/3X/0/8/082a549536bb44dd44b49da7e74906b74db37135.png)

The problem is that doing this for the 3 (and tried also with manual limits) the RGB composition does not comes very convincing. So basically, I don’t know more what to do with this type of data.