TrueColor raster satellite image

Specific Domains Geo
, ,
1

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

# 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)
2

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

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.

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

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

using GLMakie
plot(rgb)

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

1 Like
3

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?

4

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.

5

@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

6

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

:observation_data/M05

You may mean

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

7

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:

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) ;
	...
8

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.

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

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.

9

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

10

Well, I got a moderate success.

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)

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.

viz(Irgb)

GMTjl_j-or8
GMTjl_j-or81789×1806 426 KB

11

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

12

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

13

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.

1 Like
14

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
image1920×1939 737 KB

15

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?

16

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.

17

Oh right they actually nest that deeply in practice?

Damn. We probably need a new type for this.

18

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.

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

19

To be clear thats not documented API sytax :wink:

1 Like
20

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)

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
GMTjl_j1803×1856 80.3 KB

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.