Add mdt recipe

[Esteban82]

I add the recipe for the Mean Dynamic Topography.

Doubts/Issues:

1. What precision should I use? 0.01 m? The range of values is -1.47329998016 to +1.82179999352 (m I guess).
2. The master resolution is 7.5 m ( 0.125 degress or 450 sec). I think we can let the 10 m resolution.
3. Is there any typical cpt for this type of maps?

Sources: https://topex.ucsd.edu/pub/MSS_replace/ https://topex.ucsd.edu/pub/MSS_replace/faq_SWOT_MSS_rev_4.pdf

[WalterHFSmith]

Dear Federico,

Thanks for putting this in.

Yes, the range of the data is just a few meters.

I don’t know what the noise level in this data set is, but we should store it with precision of 1 mm or better, and perhaps with as much horizontal sampling as on the master. The reason I suggest these things is because the way this data is often used is as a velocity potential; that is to say, its horizontal derivatives are calculated to estimate horizontal velocity components. In the absence of noise, the true horizontal gradient of the mean dynamic topography will often be 1 micro-radian or less (1 mm per km or less). To that end, the numerically computed derivatives should not have any imprecisions contributed by an inadequate storage of the data.

I don’t have a suggestion for the typical CPT. Often a rainbow or similar is used. If you enter Mean Dynamic Topography into images.google.com http://images.google.com/ you will see many examples.

Thanks again,

Walter

On Oct 30, 2023, at 8:40 AM, Federico Esteban @.***> wrote:

I add the recipe for the Mean Dynamic Topography.

Doubts/Issues:

What precision should I use? 0.01 m? The range of values is -1.47329998016 to +1.82179999352 (m I guess). The master resolution is 7.5 m ( 0.125 degress or 450 sec). I think we can let the 10 m resolution. Is there any typical cpt for this type of maps? Sources: https://topex.ucsd.edu/pub/MSS_replace/ https://topex.ucsd.edu/pub/MSS_replace/faq_SWOT_MSS_rev_4.pdf

You can view, comment on, or merge this pull request online at:

https://github.com/GenericMappingTools/gmtserver-admin/pull/247

Commit Summary

3399f6f https://github.com/GenericMappingTools/gmtserver-admin/pull/247/commits/3399f6fe0abf037d1b94845cdb5e07eaa8d42677 add MDT recipe cab2a59 https://github.com/GenericMappingTools/gmtserver-admin/pull/247/commits/cab2a59280f54e7116ce8404af620c3d5156c659 Change year File Changes (1 file https://github.com/GenericMappingTools/gmtserver-admin/pull/247/files) A recipes/earth_mdt.recipe https://github.com/GenericMappingTools/gmtserver-admin/pull/247/files#diff-cf2df39ce41e85e5d2c64c4350fa0d1c1d533773f88d1719b8d30c28c2dd327e (38) Patch Links:

https://github.com/GenericMappingTools/gmtserver-admin/pull/247.patch https://github.com/GenericMappingTools/gmtserver-admin/pull/247.diff — Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GNZLKQAE3SMFL6KZI3YB6N2BAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43ASLTON2WKOZRHE3DQMRWGEYDSNI. You are receiving this because your review was requested.

[Esteban82]

I don’t know what the noise level in this data set is, but we should store it with precision of 1 mm or better,

I think that we have three options for the precision:

1. Use 1 mm (scale 0.001) for a 16-bit integer grid (type ns).
2. Use 0.2 mm (scale 0.0002) for a 16-bit integer grid (type ns).
3. It is not possible to use a precision 0.1 mm (or higher) for a 16-bit integer grid. For this, we should use a 32-bit float grid (nf, the same as the original). I think that it is possible becuase the original grid is not to big.

Yes, we are using the original horizontal resolution (7.5 m or 450 sec).

BTW, this is a map from the original grid.

[PaulWessel]

Since files on server are stored as JP2000 tiles we have to use integers, i.e. short int 2 byte. Given the range I think we could do scale = 1/9900 and offset = 0.17425000668 and that gives one integer equals 0.00010101010101 meter or 0.1 mm.

[PaulWessel]

So perhaps make these changes to the recipe file:

# We use a precision of 0.00010101010101  m
# The range of -1.47329998016 to +1.82179999352 m means we may use offset of 0.17425000668 and scale of 0.00010101010101
...
# DST_SCALE=0.00010101010101
# DST_OFFSET=0.17425000668

[Esteban82]

Pauil, I just made those modifications.

I am not sure if this is ok: # SRC_TITLE=SIO_Mean_Dynamic_Topography

[PaulWessel]

I think that is just fine.

[WalterHFSmith]

What you are doing sounds good to me.

On Oct 30, 2023, at 11:49 AM, Paul Wessel @.***> wrote:

I think that is just fine.

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785509455, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GICUVITG24R42R46C3YB7EBJAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVGUYDSNBVGU. You are receiving this because your review was requested.

[Esteban82]

Ok, I will process the grid to see if there is any issue.

[PaulWessel]

Thanks, yes it is a bit unusual but 1/scale is an integer (9900) but since we store the increment just make sure you stuff enough digits in the scale and offset.

[remkos]

Why not use a scale of 0.0001, so the vertical resolution is 0.1 mm? The offset can remain 0. That allows a range of ±3.2767 meters with 16-bit integer. So would easily fit.

[remkos]

What you are doing sounds good to me. … On Oct 30, 2023, at 11:49 AM, Paul Wessel @.***> wrote: I think that is just fine. — Reply to this email directly, view it on GitHub <#247 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GICUVITG24R42R46C3YB7EBJAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVGUYDSNBVGU. You are receiving this because your review was requested.

Not to me though!

[PaulWessel]

Wont fit I think. I you divide the range by 0.0001 you get 32950 which exceeds 32767:

gmt math -Q 1.82179999352 -1.47329998016 SUB 0.0001 DIV =
32950.9997368

[PaulWessel]

Would prefer 0.0001 of course but...

[remkos]

Wont fit I think. I you divide the range by 0.0001 you get 32950 which exceeds 32767:

gmt math -Q 1.82179999352 -1.47329998016 SUB 0.0001 DIV =
32950.9997368

You just proved it does fit. The dynamic range for a 16-bit integer is 65536 > 32951.

Don't forget that you have positive and negative numbers!

You would only need -14733 to +18218 from the full range of -32767 to +32767, leaving -32768 as default value. Of course, you could use +32767 for default as well.

[PaulWessel]

Radiation must have removed the signs... Yes you are right, but we could then do inc = 1/19800 = 5.0505050505e-5 and offset like before?

[WalterHFSmith]

Right. A signed 16-bit integer can represent from -32767 to +32767, and the data are smaller than +/- 3.2 m, so 1e-04 could be the scale factor.

As to the Title, I should have suggested including a version number. This information is going to evolve, and within the next year or so, because the new data are going to start pouring in.

W

On Oct 30, 2023, at 1:04 PM, Remko Scharroo @.***> wrote:

Wont fit I think. I you divide the range by 0.0001 you get 32950 which exceeds 32767:

gmt math -Q 1.82179999352 -1.47329998016 SUB 0.0001 DIV = 32950.9997368 You just proved it does fit. The dynamic range for a 16-bit integer is 65536 > 32951

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785675562, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GNYYJ5ETWLDRFQ2DKLYB7MZ3AVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVGY3TKNJWGI. You are receiving this because your review was requested.

[PaulWessel]

I will let @remkos and @walter decide. Usually we pick scale/origin to enables the highest resolution (i.e., smallest increment) as possible - hence the 0.0000505050505 vs 0.0001. This is meant to get better derivatives if you need then even though the data resolution is nowhere near 0.000050505 m.

[WalterHFSmith]

On the one hand it may be a good idea to enable the highest possible resolution. On the other hand, the maximum and minimum of the data set may change slightly as the version numbers change, and if so each version number would have to have a different scale and offset. So comparing versions numerically could also get messy. Therefore if there is a logical choice that doesn’t throw away a lot of available dynamic range, but also leaves things in a fairly convenient unit (e.g. 0.1 mm with zero offset ?) then that might be a good choice. Also, if numerical accuracy is the goal, would we have to consider the floating point representation of the scale factor? For example, would we get the most accurate results if the scale factor were a power of two or had a simple binary float representation?

W

On Oct 30, 2023, at 1:13 PM, Paul Wessel @.***> wrote:

I will let @akshmakov https://github.com/akshmakov and @walter https://github.com/walter decide. Usually we pick scale/origin to enables the highest resolution (i.e., smallest increment) as possible - hence the 0.0000505050505 vs 0.0001. This is meant to get better derivatives if you need then even though the data resolution is nowhere near 0.000050505 m.

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785690923, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GLR5WBGTSNT4HDMV3LYB7N3RAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVGY4TAOJSGM. You are receiving this because your review was requested.

[PaulWessel]

That is fine with me as a nonuser of this data set. I can offer 1/16000 = 0.0000625 meter as a cleaner increment than earlier and it should hold up for at least a few versions? But if you prefer 0.0001 just confirm that so poor @Esteban82 can finalise it.

[WalterHFSmith]

If it fits the dynamic range to use 1/16384 then that would be a better increment, no?

Federico, please use 1/16384 if it fits, otherwise use 1/10000. Thanks, W

On Oct 30, 2023, at 2:28 PM, Paul Wessel @.***> wrote:

That is fine with me as a nonuser of this data set. I can offer 1/16000 = 0.0000625 meter as a cleaner increment than earlier and it should hold up for at least a few versions? But if you prefer 0.0001 just confirm that so poor @Esteban82 https://github.com/Esteban82 can finalise it.

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785811817, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GIYOWP3HVCTRKJG7WLYB7WUHAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVHAYTCOBRG4. You are receiving this because your review was requested.

[PaulWessel]

OK, but 1/16384 = 6.103515625e isn't cleaner (fewer decimals needed) than 1/16000 = 6.25e-05 but it undoubtedly smaller :-).

[WalterHFSmith]

I’m not picky on this. But I thought 1/16384 could be represented exactly in binary floating point, whereas 1/16000 or 1/10000 cannot. Is that correct?

W

On Oct 30, 2023, at 2:45 PM, Paul Wessel @.***> wrote:

OK, but 1/16384 = 6.103515625e isn't cleaner (fewer decimals needed) than 1/16000 = 6.25e-05 but it undoubtedly smaller :-).

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785838867, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GMLIPDLKRIS2K6FRXTYB7YVPAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVHAZTQOBWG4. You are receiving this because your review was requested.

[PaulWessel]

I think that is correct. It just means we must do lots of decimals in

# DST_SCALE=0.00006103515625....

to approximate 1/16384 but do we really get that?

[WalterHFSmith]

I don’t know how this is used inside our machinery. If we have to specify it as some character string that is going to get converted to floating point, then that won’t work. But if we can specify something like

z_to_be_stored = (short) rint(z * 16384); z_scale_factor = 1 / 16384.0;

Then we would be ok.

If this isn’t possible, and we are writing ascii look-up tables with character strings that are scanned to floating point numbers, then might has well just use 1/10000

W

On Oct 30, 2023, at 2:53 PM, Paul Wessel @.***> wrote:

I think that is correct. It just means we must do lots of decimals in

DST_SCALE=0.00006103515625....

to approximate 1/16384 but do we really get that?

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785850938, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GNO7PXL3ZUJMULFO5TYB7ZT5AVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVHA2TAOJTHA. You are receiving this because your review was requested.

[PaulWessel]

Yes, we have an ASCII table that every user downloads automatically and it is ASCII. There is no netCDF file until we convert the jp2 tiles to netCDF relaying on the meta data in the ASCII table. We also have bash scripts that does all the processing and filtering and downsampling. The objective is to make the tiles as small as possible for global download. We do this by using JP2 since that is the max compression we can get.

[remkos]

I would definitely stick with 0.1 mm. I haven't looked into the original file, but if it is any similar to the previous CNES/CLS solution (MDT CNES/CLS 2022), then the original grid already had a vertical resolution of 0.1 mm ( in that case needlessly using a 32-bit int)

[Esteban82]

Since files on server are stored as JP2000 tiles we have to use integers, i.e. short int 2 byte.

Just a comment. The highest resolution is 7.5 m. Therefore, none of the grids will be tiled. At least for this version.

So,. we could use a 32-bit float.

[WalterHFSmith]

So, if no tiling then question of scaling is moot for now?

w

On Oct 30, 2023, at 3:27 PM, Federico Esteban @.***> wrote:

Since files on server are stored as JP2000 tiles we have to use integers, i.e. short int 2 byte.

Just a comment. The highest resolution is 7.5 m. Therefore, none of the grids will be tiled. At least for this version.

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/gmtserver-admin/pull/247#issuecomment-1785899876, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GO2LGJ5WN3R4P3DCPDYB75RZAVCNFSM6AAAAAA6V6JYM6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBVHA4TSOBXGY. You are receiving this because your review was requested.

[PaulWessel]

Heh, I somehow missed the resolution. Yes, only 5m and higher resolutions are tiled, so for now no tiles. Suspecting you will push that resolution up as more data comes along and then we can revisit this chat!

[Esteban82]

Ok, so I have use this parameters

# DST_FORMAT=ns
# DST_SCALE=0.0001
# DST_OFFSET=0

I think that the original source is: Jousset S., Mulet S., Wilkin J., Greiner E., Dibarboure G. and Picot N.: “New global Mean Dynamic Topography CNES-CLS-22 combining drifters, hydrological profiles and High Frequency radar data”, OSTST 2022, https://doi.org/10.24400/527896/a03-2022.3292 .

[Esteban82]

Paul, when I run srv_tiler.sh I got that message.

bash scripts/srv_tiler.sh earth_mdt
No tiling requested for earth_mdt_07m_g.grd
scripts/srv_tiler.sh: line 220: printf: 7.5: invalid number
scripts/srv_tiler.sh: line 220: printf: 7.5: invalid number
No such file to tile: /home/federico/.gmt/server/earth_mdt_07m_p.grd

In the earth_mdt_server.txt error I got this info (it says 07m instead of 7.5m).

/server/earth/earth_mdt/ earth_mdt_07m_g.grd 07m g 0.0001 0 2.8M 0 2023-10-30 - - @earth_mdt.cpt SIO Mean Dynamic Topography at 7x7 arc minutes reduced by Gaussian Cartesian filtering (39.3 km fullwidth) [Sandwell et al., 2023]

[PaulWessel]

Hm, yes 7.5m is not going to work unless we update some scripts, since naming is planet_dataset_##m_g|p.grd . Give this a try. It places 7.5m in the meta file and hopefully that is all we need. The 7.5m file will be called 07m but it will be 7.5m as far as increments go. I will look at this more after my doctor apt.

[Esteban82]

I change the title to Mean_Dynamic_Topography_CNES_CLS22. I think it is more precise (and also includes the version).

[Esteban82]

Give this a try.

I try it. In the earth_mdt_server.txt, there info of the grid earth_mdt_07m_p.grdwhich doesn't exists (only the _g version is).

/server/earth/earth_mdt/ earth_mdt_07m_p.grd 7.5m p 0.0001 0 2.8M 0 2023-10-31 - - ...

[PaulWessel]

Deep in pyGMT debugging right now so need to finish that. The script should check if master and then it will not try to make the other resolution. Perhaps confusion with 07m file name? Run the script with bash -xv to see what is happening in the loop.

[PaulWessel]

I think it is not worth trying to code up a scheme that skips outputting the nonexistent 07m_p in the txt file. They will soon enough up the resolution with the next version (data is coming in as we speak) so let us just remember this, remove the offending line in earth_mdt_server.txt, and upload to candidate. Just done that so mdt is on candidate and no 07m_p in the table.

[PaulWessel]

I pushed some changes to Makefile so it knows about the new data set.

[PaulWessel]

BTW, getting lots of this when the final chmod command in make place-earth-mdt, e.g.

chmod: changing permissions of '/export/gmtserver/gmt/candidate/server/earth/earth_mss/earth_mss_01m_g/S60W150.earth_mss_01m_g.jp2': Operation not permitted

Perhaps your (and my) umask prevents the other from running chmod?

[PaulWessel]

With no cpt, do we create one based on turbo and add it to the cache as earth_mdt.cpt?

[Esteban82]

I think that turbo is ok. An alternative would be a CPT like polar with a hinge in 0. But, I think turbo is better.

[PaulWessel]

Would you mind making a suitably ranged earth_mdt.cpt, add it to gmtserver-admin/cache and do a PR. Takes at least 1 hour for the server to sync that up once we merge. Then maybe add mdt to series of plots?

[Esteban82]

I think it is not worth trying to code up a scheme that skips outputting the nonexistent 07m_p in the txt file. They will soon enough up the resolution with the next version (data is coming in as we speak) so let us just remember this, remove the offending line in earth_mdt_server.txt, and upload to candidate. Just done that so mdt is on candidate and no 07m_p in the table.

I saw that you upload the files. I think that the problem is that the script (srv_tiler.sh) can't see the record that has the master code.

# res   unit    tile    chunk   code
03      m       90      2048    master
04      m       180     2048

See that in the file earth_mdt_server.txt says: `earth_mdt_07m_g.grd ... at 7.5x7.5 arc minutes reduced by Gaussian Cartesian filtering (39.3 km fullwidth)' which is not true for the original resolution.

[Esteban82]

I found out that the script don't see the master code if there is no tiling.

I ask for a tile of 360 in the recipe and I got the correct records

# List of desired output resolution and chunk size.  Flag the source resolution with code == master
# res   unit    tile    chunk   code
7.5     m       360     4096    master

Record for earth_mdt_server.txt:

/server/earth/earth_mdt/ earth_mdt_7.5m_g/ 7.5m g 0.0001 0 2.8M 360 2023-11-03 - - @earth_mdt.cpt Mean Dynamic Topography CNES CLS22 original at 7.5x7.5 arc minutes [Jousset et al., 2022]

[PaulWessel]

OK, maybe edit the *txt to say things similarly to other master datasets elsewhere?

[PaulWessel]

So that gives you a tile directory with 1 tile, then? I guess that is fine if that improves the metadata.

[Esteban82]

So that gives you a tile directory with 1 tile, then? I guess that is fine if that improves the metadata.

With 360 I got an empty directory. The original file was moved to earth_mdt_tiled.

[PaulWessel]

OK, I think I can modify the script so that if tile is 360 is a special message that says do not move the file to the tiled directory but treat as other single files?

[Esteban82]

I think that can work. Although it would be strange if it said 360 in the recipe.

Another option would be for the script not to pay attention to the tiling when writing the information about the master grid.

[Esteban82]

BTW, I think that this PR can be approve.