Harpmerge on the command line | ingestion question on GOME2AB SO2 flags

[mariliza]

Dear Sander, dear all forum users,

I wish to make some daily gridded files for the GOME2AB SO2 operational GDP4.8 product. I am using :

harpmerge -o ‘detailed_results=SO2;corrected=true;so2_column=2.5km’ -ap ‘bin_spatial(361,-90.,0.5,721,-180.,0.5);derive(latitude {latitude});derive(longitude {longitude});derive(SO2_column_number_d
ensity [DU]);derive(SO2_slant_column_number_density [DU]);exclude(weight,solar_zenith_angle_weight,SO2_column_number_density_count)’ -a ‘keep(latitude_bounds,longitude_bounds,SO2_column_number_dens
ity,solar_zenith_angle, cloud_fraction,SO2_column_number_density_validity, scan_direction_type, SO2_slant_column_number_density); solar_zenith_angle <= 75; cloud_fraction <= 0.2; scan_direction_ty
pe==0’ -l $filename outputfile.nc

Is there a way to filter only for the SO2_flag using harpmerge as such, via the command line and not via python i.e. using both harp ingest and coda_fetch and then merging? for this product, for the 2.5km placing SO2_column_number_density_validity==0 is fine as I am not interested in the volcanic component. But when I use so2_column=15km I also wish the volcanic pixels left untouched, which requires a different kind of filtering.

Many thanks!
MariLiza

[sander.niemeijer]

I am not entirely sure what you mean, do you mean just adding a SO2_column_number_density_validity==0 filter?

If so, you can just add it to the start of the -a operation block:

harpmerge -o ‘detailed_results=SO2;corrected=true;so2_column=2.5km’ -ap ‘bin_spatial(361,-90.,0.5,721,-180.,0.5);derive(latitude {latitude});derive(longitude {longitude});derive(SO2_column_number_density [DU]);derive(SO2_slant_column_number_density [DU]);exclude(weight,solar_zenith_angle_weight,SO2_column_number_density_count)’ -a ‘SO2_column_number_density_validity == 0; keep(latitude_bounds,longitude_bounds,SO2_column_number_density,solar_zenith_angle, cloud_fraction,SO2_column_number_density_validity, scan_direction_type,SO2_slant_column_number_density); solar_zenith_angle <= 75; cloud_fraction <= 0.2; scan_direction_type==0’ -l $filename outputfile.nc

[mariliza]

Apologies. What I mean is this: according to the PUM, Microsoft Word - DLR_GOME-2_PUM_3B_Nov2019 (acsaf.org), on page 46, it is stated that there exists an SO2_flag [Table 6.7.6] whose value, if set to 0, assures no SSA, no polar regions and other nastiness. There exists also an SO2_Volcano_Flag [Table 6.7.7] which assumes a value of 0 when no volcano is detected.

Now, if I want to access the SO2=15km product, i.e. the volcanic product, if I set SO2_column_number_density_validity==0 I also exclude the volcanic flagged pixels, see the two [quick] plots attached. In the one without filtering, the Hawaii emissions are there, in the one with SO2_column_number_density_validity==0 the Hawaii eruption disappears.

Hope I explained it better!
M

ariLiza

[sander.niemeijer]

Ok. I think I start to understand what you mean. The SO2_column_number_density_validity variable in HARP is a combination of the three flag variables QualityFlags, SO2_Flag, and SO2_Volcano_Flag from the product. These are treated as a bitfield.

If you only want to have values where SO2_Flag==0 (and QualityFlags and SO2_Volcano_Flag can be any value), you would have to use a bitfield operation on it. Since SO2_Flag is bit 5-8 within SO2_column_number_density_validity, you would have to match against 2^4 + … + 2^7 = 240.
So, if you use SO2_column_number_density_validity !& 240 you will only get measurements for which bits 5-8 are 0.

I believe this is what you were looking for.

[mariliza]

Good morning Sander,

I think we are now getting somewhere!

You load onto “SO2_column_number_density_validity” the Quality Flag, as bits 0-3, the SO2_Flag, as bits 4-7 and the SO2_Volcano_Flag as bits 8-11. correct?

For e.g. if I want to be sure that the retrieval is overall fine and it is a volcanic eruption, I need bits 0-7 not raised and bits 8 & 9 to be raised, i.e. 2^8 + 2^9 = 768 will be my check?

Thanks again,
MariLiza

[sander.niemeijer]

@mariliza

You load onto “SO2_column_number_density_validity” the Quality Flag, as bits 0-3, the SO2_Flag, as bits 4-7 and the SO2_Volcano_Flag as bits 8-11. correct?

Almost. See the mapping documentation. QualityFlags is bits 0-3, SO2_Flag is bits 4-7, but SO2_Volcano_Flag is mapped using 2^(SO2_Volcano_Flag-1) which only captures 3 bits (bits 8-10).

@mariliza

For e.g. if I want to be sure that the retrieval is overall fine and it is a volcanic eruption, I need bits 0-7 not raised and bits 8 & 9 to be raised, i.e. 2^8 + 2^9 = 768 will be my check?

If you use SO2_column_number_density_validity =& 768 this will return true if only bits 8 and 9 are set, but it will not care about the values of bits 0…7. If you want bits 0…7 to be 0 and bits 8 and 9 to be 1, you will have to perform two separate filters: SO2_column_number_density_validity !& 255; SO2_column_number_density_validity =& 768.

[mariliza]

Dear Sander,

I am honoured that you think I can simply read and understand this,

out of its computational context!

Now that you’ve explained this in English, yes, I follow. Careful thinking is required to make sure I do not filter too much.

How about adding this thread as a hyperlink to the mapping documentation?! I am sure it will help more than one person in properly filtering their SO2 fields.

Best wishes
MariLiza

[mariliza]

Dear Sander,

I have now come back to this issue and have a further clarification that I would greatly appreciate. I start with what I have come to understand from our chat and by reading the PUM and by reading the data both using harp as well as the simple netCDF4 python module.

The original SO2_Volcano_Flag is not in bits, but just one number which may be either 0, 1, 2 or 3, as per the PUM [https://acsaf.org/docs/pum/Product_User_Manual_NTO_OTO_Nov_2019.pdf].

The SO2_column_number_density_validity` variable in HARP is a bitfield, and you have translated those 4 possible values into three bits, namely bit 8, 9 and 10, correct? so you are ignoring the value of 0, obviously, since it means no volcanoes, and giving the opportunity for someone to check if the values of 1, 2 or 3 appear.

Now, contrary to bits 0 - 7, which we want to be zero, i.e. not raised, for a volcanic pixel we require either bit 8 [original value 1] OR bit 9 [original value 2] to be raised. Not both at the same time, either/or.

I went through the bitfield operators HAPR has, but I fail to see this either/or option. Might have been an option on numerical values [<=, >=] but not on bitfields. If I understand the online harp explanations on the bitfield operators properly.

Any way out of this predicament?

Of course I can always simply ignore bit 9, original value 2, which tells me there is volcano near an anthropogenic source. If the PUM speaks the truth then all original value 2 SO2_Volcano_Flag pixels have also the original value 1 SO2_Volcano_Flag pixels included. I will not be missing anything.

Many thanks, as ever,
MariLiza

[sander.niemeijer]

I should have realised this before. You can never have both bits 8 and 9 set at the same time due to the way we map the flag.

We didn’t actually have an operator in HARP to filter for ‘any’ bits in the mask to be non-zero.
I just added this to HARP (the ‘=|’ operator). This will be included in the next HARP release.

Just for me to test this a bit. Do you know of a GOME2 orbit where there is actually an SO2_Volcano_Flag set to 2? I checked for instance 21 Sep 2021, where S5P shows a plume going from Sicily to Greece, where I would potentially expect such a flag near Greece, but GOME2(b) doesn’t show anything in that area (maybe due to different overpass time?)

[mariliza]

Hi Sander,

Yes, you can try GOME2A 2010 05 15, it was during the Eyja eruption. In one of the orbits I found a volc flag of 2, and in multiple others volc flag of 1. Should I try and find out in which ones?

Just so you can verify, this one is with SO2_column_number_density_validity !& 255 only:

And this one, also with SO2_column_number_density_validity =& 256

Thanks!

[sander.niemeijer]

I got all the O3MOTO products from acsaf.eoc.dlr.de for GOME-2A for 2010-05-15, but I only found SO2_Volcano_Flag values 0 and 1.

[mariliza]

Oh, of course! I am working on the reprocessed dataset. Check your email.

[sander.niemeijer]

Thanks for the test data. I was able to verify that using SO2_column_number_density_validity=|768 now does the trick with new HARP change.

[mariliza]

Thanks Sander. Will you include it in v1.17? Release 1.17 · stcorp/harp · GitHub this is also the one avaliable via conda install -c conda-forge harp.

I note that it exists as one of the two commits here: Comparing 1.17…master · stcorp/harp · GitHub

But it is not clear to me if these are also in the anaconda channel?

Thanks!
MariLiza

[sander.niemeijer]

A version of HARP that is released will not get modified (that is fundamentally how versioning works).

As I mentioned, this new operator will be included in the next HARP release. When that will arrive exactly isn’t planned yet.

Users that really need a recent change that was committed can install HARP from github (which is a bit more involved than just performing a conda install).

[mariliza]

Hey Sander,

I do not wish to be the bearer of bad news, however I should probably warn you that there is a difference between the volcanic flag for GOME2A & 2B, compared to GOME2C. If you look at page 46-47 of the PUM, Microsoft Word - DLR_GOME-2_PUM_3B_Nov2019 (acsaf.org), you will note that there is an extra value for the GOME2C products.

Not sure if this eventuality is covered by the new commit you made last month.

Many thanks,
MariLiza

[sander.niemeijer]

We were not aware of this. It seems rather inconvenient that they chose to use different processor versions for the different MetOp satellites. This makes it rather hard to combine MetOp B and C SO2 data.

I will see whether we can accommodate these changes.

[sander.niemeijer]

I checked, and I think we can leave the current HARP implementation as-is. It will then be up to the user to interpret the resulting SO2_column_number_density_validity variable appropriately according to whether it comes from MetOp B or C (as they would have had to do if they would have extracted the SO2_Volcano_Flag directly from the product instead of using HARP).

I don’t see an easy way to ‘harmonise’ the flags between GDP 4.8 and 4.9.

[mariliza]

Thanks for checking, Sander! Indeed, it makes using the command line harp tools more complicated. I should mention here that yours truly [and team!] is currently validating the reprocessed GOME2A & 2B which are also run with GDP4.8, and not GDP4.9, so the issue will be with us for some time to come.

On the positive side, the extra volcanic flag of GDP4.9 can be tackled with in GDP4.8 using the flag that covers the SAA region with NaNs, and of course adding a SZA filter is always a good idea with SO2.

Many thanks for all,
MariLiza