Lat and lon coordinates empty in output file

[davidorme]

The output netCDF file does not have the lat and lon coordinates set correctly. Using xarray in python to explore the data, I get this output:

n [87]: out                                                                                                                                 
Out[87]: 
<xarray.Dataset>
Dimensions:  (lon: 720, lat: 360, pft: 13)
Coordinates:
  * lon      (lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
  * lat      (lat) float32 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
  * pft      (pft) int32 1 2 3 4 5 6 7 8 9 10 11 12 13
Data variables:
    biome    (lat, lon) float32 27.0 27.0 27.0 27.0 27.0 ... 0.0 0.0 0.0 0.0 0.0
    wdom     (lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    gdom     (lat, lon) float32 279.0 279.0 279.0 279.0 ... 0.0 0.0 0.0 0.0
    npp      (pft, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 ... nan nan nan nan nan
Attributes:
    title:        BIOME4 netCDF output file
    timestamp:    20220707 1205
    Conventions:  COARDS
    node_offset:  1

In [88]: for v in out.variables:  
    ...:     print(f"{v}: {np.nanmin(out[v])} to {np.nanmax(out[v])}")  
    ...:                                                                                                                                     
lon: 0.0 to 0.0
lat: 0.0 to 0.0
pft: 1 to 13
biome: 0.0 to 27.0
wdom: 0.0 to 14.0
gdom: 0.0 to 9999.0
npp: -6.561651222369846e-11 to 3679.67529296875

Actually - come to think of it gdom looks odd too!

[jedokaplan]

Can you try a run using the "standard" input file I posted on the other issue thread and report back?

[davidorme]

I think the issue is that it the lat and lon values are not read from the input file, but are allocated:

https://github.com/jedokaplan/BIOME4/blob/ff3d4ff2d64568d717d71c0d295fc96400c09dc0/biome4driver.f90#L124-L125

Presumably, that allocation either has a default 0 value or we've just seen it being allocated 'empty' memory.

At the moment, I think the coordstring subsetting expects array indices, rather than allowing GMT strings with lat and long values ( e.g. -R-180/180/-90/90). That would obviously require more logic to identify the index on lat and lon dimensions.

[davidorme]

I think I've fixed this up:

• I've used the existing functions calcpixels and nearest in coordsmod.f90.
• I have renamed nearest to nearest_coords to avoid shadowing an existing name and also added include_coords as an alternative to ensure that provided GMT coordinates are included (where possible) within the the subset.
• The source lat and lon data is now read into two new vectors (all_lat and all_lon) and then the appropriate slice is extracted into the lat and lon vectors, given the subset indices.