sakov / enkf-c Goto Github PK

Hi Pavel,

I'm wondering about the writing/reading/writing scheme in enkf_prep.c.

Why prep has to re-read the superobservations file and discard the sobs structure that was just used to write the superobs file?
Also, why the obs_checkforland has to be run after writing and not before?

Finally, shouldn't the write of sobs file be postponed to be after land checking and avoid write/read/write thing? I did some changes to "write" at the end only, and the observations.nc are unchanged for my tests. I'm wondering if this write/read/write is some "cleanup" on the sobs structure or just a remnant of testing...

    if (do_superob) {
        enkf_printf("  superobing:\n");
        obs_superob(obs, cmp_obs, &sobs, describe_superob_id);

        if (describe_superob_id >= 0)
            goto finalise;

        enkf_printf("  writing superobservations to \"%s\":\n", FNAME_SOBS);
        obs_write(sobs, FNAME_SOBS);
        // WHY!?
        free(sobs->data);
        enkf_printf("  reading super-observations from disk:\n");
        obs_read(sobs, FNAME_SOBS);

        enkf_printf("  checking for superobs on land:\n");
        if (obs_checkforland(sobs, m)) {
            obs_compact(sobs);
            for (i = 0; i < sobs->nobs; ++i)
                if (sobs->data[i].status != STATUS_OK)
                    break;
            assert(i != sobs->nobs);
            enkf_printf("    deleted %d observations\n", sobs->nobs - i);
            sobs->nobs = i;
            //WHY bother!?
            file_delete(FNAME_SOBS);
            enkf_printf("  writing good superobservations to \"%s\":\n", FNAME_SOBS);
            //move outside the if?
            obs_write(sobs, FNAME_SOBS);
        } else
            enkf_printf("    all good\n");
    } else {
        observation* data = malloc(obs->ngood * sizeof(observation));

        memcpy(data, obs->data, obs->ngood * sizeof(observation));
        sobs = obs_create_fromdata(obs, obs->ngood, data);
        obs_write(sobs, FNAME_SOBS);
        goto finalise;
    }

Hi pavel,

I would like your opinion about some testing I did regarding the bin averaging from enkf_prep.

I'm having two behaviors that are related and unexpected at first sight. They seem to be a small bug or an approximation, but for me is hard to tell today if any adjustment is done within the enkf-c kernel to account for the following:

1. The superobs is ignoring the boundary points. This is very clear in the following figure:

The above is the lower/left boundary of the "EAC" domain (example 1). The model boundary is mark with a blue line, The raw observations (blue) is a high res. sst gridded data, the superobs arr in red, while the center of the cell (grid-t) is a star, with the cell-edges in gray. (see link at the bottom to obtain the folder required to run the example).

Clearly the boundaries are being ignored. If one want to update the boundary estimates, the code will not use local obs (in the grid sense), only through localization and nearby points.

This plot also show the second problem:

2. The superobs locations are shifted in physical space by half a grid-spacing. Thus, they are not where they supposed to be in regard to the grid.prm coordinates and state variable definitions.

I would expect that, if I want superobs of SST, they should be bin averaged to the center of the cell (t-grid), not the corners (as indicated above). However, ENKF-C put the superobs at the corner of the cells not in their center.
I'm assuming the code is using the dx/dy on grid-t (center) to estimate the cell width/height but it did not account for the first points and do not adjust it's position (half-shift to the left and below). This is clear by the lack of boundary points (as above) and the aproximate locations of superobs at the actual corners:

On the above, this behaviour clearly generate wrong positions of superobs. In fact, the bin averaging are shifiting the "true" field up and to the right -- because it's assuming the cell corners are the actual centers. This become even more problematic in curvilinear grids:

On the figure above, it's clear that ENKF-C is doing an average where is not required -- there are superobs between points that are actually inside a single cell and very close to the actual center.

Clearly this is a matter of grid stagering location and dx/dy of the cells. If I provide grid.prm with the corner locations (psi points in ROMS for example), I got the expected results:

Again, using the corners as input to enkf_prep I got correct positions, since now the dx/dy of these positions are at the actual center of the cell. (I still got no points at the boundary but this is because the reduced index space of the cell corner array).

So this is expected? Should be disregarded/assumed a good aproximation given all the examples in enkf-c provide center cells but superobs are not actually there!?

At first sight appears only a slight correction is needed (shift the positions/dx/dy to the left/down in the code).

I can upload the prm/nc files required if this is indeed a bug.

The current generic readers do not support several flags to be used as valid data in obs.prm

second commit in #11 add/fix support to:

parameter qcflags = 1, 2, 3, 4

to consider the flags as good ones.

The current logic on estd for the generic reader is incomplete:

1. If the user did not provide an estd variable (parameter estdname is missing)
2. but provide an std variable (parameter stdname = VARSTD)
3. The code completely miss setting o->std with the proper values in stdname. EDITED

See line 361 onwards in reader_xy_gridded.c for example.

The third commit in #11 fix this issue by checking if both std/estd where not provided, before setting the std and keep going.
The check inside is pedantic on purpose ( the nan check), given that I don't know if you want to set a precedence order in the case the user provided all three std parameters (estd/std/ and/or error_std attribute).

The last commit in #11 fix a segfault case where estd[i] is accessed when not defined because the change in logic above.

Hi,

I have been trying to compile enkf library, but i keep on getting this error

common/observations.c:878:37: error: unknown type name '__compar_d_fn_t'
void obs_superob(observations* obs, __compar_d_fn_t cmp_obs, observations** sobs, int sobid)

I think it is an incompatibility issue between mac and linux. I can compile the rest of files by commenting out "__compar_d_fn_t" usage. Also, when i compile, the qsort-r gets called from stdlib, but i think i still need to have way around this as this is used by super-orb part of this library

i would appreciate any help with this.

for running make enkf, I am getting the following error.
error is:
writing point logs:
calculating transforms for log point (156.000,-32.000): 539 obs
*** buffer overflow detected ***: ./enkf_calc terminated.

I gave NP_CALC = 1
NP_UPDATE = 1
for NP_CALC=4 also getting the same error.

while compiling, I am getting the following warning
common/utils.c: In function ‘get_memory_usage_kb’:
common/utils.c:1324:12: warning: ignoring return value of ‘fread’, declared with attribute warn_unused_result [-Wunused-result]
(void) fread(buffer, sizeof(char), to_read, procfile);

Hi,

I tried to execute all provided examples with provided datasets which ran quite easily. But i noticed that .nc files are just holding metadata which is appearing in Panoply (NC reader) tool but nothing to plot or see in these files. Can you please share the reason of this strange type of files. secondly, How can I see the results as there are number of files are generated by examples script under their respective folders?

I also tried to go through the Guide book but that is more about theory instead of your code description. Can you please let me know if you have any detailed document about code description.

Thanks.

Hi,

I think It would be good to be able to specify all input files as arguments:

enkf_prep -i enoi.prm -obstype obstype.prm -obs obs.prm -grid grid.prm

This would enable code tests to be easily implemented and executed, such as in the makefile. This would also help identify bugs and test invalid user requests.

This led me to think on why not put almost everything in one single file!? maybe I'm missing something, but except for the obs.prm, I cannot foresee advantage of several input files. It would help a lot for testing too by simplifying the call and case writing.

Moreover, in a single file, changing parameters would be 1. easy with cat/sed/awk/etc, 2.help to maintain a consistent experiment when copying around, 3. running different parameter cases would be easier, and 4. organization would be better.

In any case, people out there are using generators of the prm files. This would only simplify their life, by including a single generator of one (bigger) input file (as well as the obs.prm).

Hi @sakov,

The function dir_rmifexists fails if the directory is not empty. The use of rmdir in https://github.com/sakov/enkf-c/blob/master/enkf/common/utils.c#L619 only works if it's an empty directory.

I have seen this impact one of our jobs that got killed in a non-clean manner such that the DIRNAME_TMP is not cleaned up. I have pasted the output below from an attempt to rerun after this event.

I think you need something like https://stackoverflow.com/questions/5467725/how-to-delete-a-directory-and-its-contents-in-posix-c/5467788#5467788 to handle recursive deletion instead.

  ERROR: enkf: CPU #0: could not remove directory ".enkftmp": Directory not empty


  I am CPU #0, now printing the backtrace:

  obtained 5 stack frames:
[0x45173d]
[0x44e824]
[0x43315a]
[0x325bc51]
[0x4081a1]

Hi Devs,

I am trying to add in the SST dataset (GHRSST) for assimilation from CMEMS database, example:
https://www.ncei.noaa.gov/data/oceans/ghrsst/L3C/GLOB/AVHRR_SST_METOP_B_GLB/OSISAF/2022/013/

I have specified in obs.prm as below:

product == GHRSST
type = SST
reader = gridded_xy
parameter VARNAME = sea_surface_temperature
parameter VARSHIFT = -273.15
file = obs/20180102120000*.nc

The program seems to be able to read it, but I am getting the following outputs

PRODUCT = GHRSST, TYPE = SST, reader = gridded_xy
reading obs/20180102120000-OSISAF-L3C_GHRSST-SSTsubskin-AVHRR_SST_METOP_B_GLB-sstglb_metop01_20180102_120000-v02.0-fv01.0.nc:
VARNAME = sea_surface_temperature
LONNAME = lon
LATNAME = lat
(ni, nj) = (7200, 3600)
nobs = 8153894
# obs added = 40276
id = 2504 - 42779
total 40276 observations
40161 valid observations
min day = 0.000
max day = 0.000
bad observation detected: type = SST, product = GHRSST, instrument = METOP_B, datafile = obs/20180102120000-OSISAF-L3C_GHRSST-SSTsubskin-AVHRR_SST_METOP_B_GLB-sstglb_metop01_20180102_120000-v02.0-fv01.0.nc, id = 2504, original id = 2504, batch = 0, value = 23.2, estd = nan, footprint = 0, lon = 96.725, lat = -2.325, depth = 0.0, model_depth = 4730.2, fi0 = 85.165, fi1 = 0.522, fi2 = 0.000, fk = 0.000, time = 0, status = 0

ERROR: enkf: CPU #0: bad observation

Have I specified the prm correctly?

Hi pavel,

Apparently, the current HEAD code is superobing at the wrong positions in the vertical. This appears similar to #5.

Looks like cell centres for sigma type coordinates are assumed as cell-edges:

The above show clearly that enkf_prep is superobing something that should not be averaged, given obs lie within the cell edges (black lines). Looks like an offset like in the one done in #5 is required. Not sure about other coordinates, but likely the case too.

Observations are from a fake z profile using the cars readers. All depths were calculated with the roms equations for Vt=2,Vs=4,theta_s=7,theta_b=0.5,hc=200. Local depth is 4883.527.

I will tear down an example if required.

Hi,

The generic readers are hardcoded to block qcflags larger than 31. I don't think this is generic enough. Several satellite products have QC flags much larger than 31 (sparse like 1,2,4,6,31,256,2048) while others are more compact (1..10).

I believe the only limiting factor here is to make sure they are integers.

See/discuss first commit in #11.