Due to the lapse in federal government funding, NASA is not updating this website. We sincerely regret this inconvenience.
A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
sssref.c
Go to the documentation of this file.
1 #include <sys/types.h>
2 #include <unistd.h>
3 
4 /* ============================================================================ */
5 /* module sssref.c - retrieve salinity from WOA climatology */
6 /* ============================================================================ */
7 
8 #include "l12_proto.h"
9 #include "h5io.h"
10 
11 static float sssbad = BAD_FLT;
12 static int32_t format = -1;
13 
14 
15 /* ----------------------------------------------------------------------------------- */
16 /* get_woasssclim() - read and interpolate WOA salinity climatology */
17 /* ----------------------------------------------------------------------------------- */
18 
19 #define WOASSSNX 360
20 #define WOASSSNY 180
21 
22 float get_woasssclim(char *sssfile, float lon, float lat, int day) {
23  static int firstCall = 1;
24  static int nx = WOASSSNX;
25  static int ny = WOASSSNY;
26  static float dx = 360.0 / WOASSSNX;
27  static float dy = 180.0 / WOASSSNY;
28 
29  typedef float sssref_t[WOASSSNX + 2];
30  static sssref_t* sssref;
31  static sssref_t* sss_sav;
32  //static float sssref[WOASSSNY+2][WOASSSNX+2];
33  //static float sss_sav[WOASSSNY+2][WOASSSNX+2];
34 
35  float sss = sssbad;
36  int i, j, ii;
37  float xx, yy;
38  float t, u;
39 
40  if (firstCall) {
41  sssref = (sssref_t*) malloc((WOASSSNY + 2) * sizeof (sssref_t));
42  sss_sav = (sssref_t*) malloc((WOASSSNY + 2) * sizeof (sssref_t));
43 
44  float ssstmp[WOASSSNY][WOASSSNX];
45  char name [H4_MAX_NC_NAME] = "";
46  char sdsname[H4_MAX_NC_NAME] = "";
47  int32 sd_id;
48  int32 sds_id;
49  int32 rank;
50  int32 nt;
51  int32 dims[H4_MAX_VAR_DIMS];
52  int32 nattrs;
53  int32 start[4] = {0, 0, 0, 0};
54  int32 status, ix, iy, ct;
55  int32 lymin, lxmin, lymax, lxmax;
56  float sum;
57 
58  int16 mon = (int) day / 31 + 1; // month of year (no need for perfection)
59 
60  firstCall = 0;
61 
62  printf("Loading SSS reference from Climatology file: %s\n", sssfile);
63  printf("\n");
64 
65  /* Open the file */
66  sd_id = SDstart(sssfile, DFACC_RDONLY);
67  if (sd_id == FAIL) {
68  fprintf(stderr, "-E- %s line %d: error reading %s.\n",
69  __FILE__, __LINE__, sssfile);
70  exit(1);
71  }
72 
73  /* Get the SDS index */
74  sprintf(sdsname, "sss%2.2i", mon);
75  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
76  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
77  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) & ssstmp[0][0]);
78  if (status != 0) {
79  fprintf(stderr, "-E- %s Line %d: Error reading SDS %s from %s.\n",
80  __FILE__, __LINE__, sdsname, sssfile);
81  exit(1);
82  }
83  status = SDendaccess(sds_id);
84  status = SDend(sd_id);
85 
86  /* rotate 180-deg and add wrapping border to simplify interpolation */
87  /* new grid is -180.5,180.5 in i=0,361 and -90.5,90.5 in j=0,181 */
88 
89  for (j = 0; j < ny; j++) {
90  for (i = 0; i < nx; i++) {
91  ii = (i < nx / 2) ? i + nx / 2 : i - nx / 2;
92  if (ssstmp[j][i] > 0)
93  sssref[j + 1][ii + 1] = ssstmp[j][i];
94  else
95  sssref[j + 1][ii + 1] = sssbad;
96  }
97  sssref[j + 1][0] = sssref[j + 1][nx];
98  sssref[j + 1][nx + 1] = sssref[j + 1][1];
99  }
100  for (i = 0; i < nx + 2; i++) {
101  sssref[0] [i] = sssref[1][i];
102  sssref[ny + 1][i] = sssref[ny][i];
103  }
104  /*
105  * Extend the grid by 1 point (use only full grid boxes later)
106  */
107  memcpy(sss_sav, sssref,
108  (WOASSSNY + 2) * (WOASSSNX + 2) * sizeof ( float));
109  for (iy = 0; iy < ny + 2; iy++) {
110  lymin = (iy == 0) ? 0 : iy - 1;
111  lymax = (iy == ny + 1) ? ny + 1 : iy + 1;
112 
113  for (ix = 0; ix < nx + 2; ix++) {
114  if (sssref[iy][ix] < sssbad + 1) {
115  lxmin = (ix == 0) ? 0 : ix - 1;
116  lxmax = (ix == nx + 1) ? nx + 1 : ix + 1;
117 
118  sum = 0.;
119  ct = 0;
120  for (j = lymin; j < lymax + 1; j++)
121  for (i = lxmin; i < lxmax + 1; i++) {
122  if (sss_sav[j][i] > sssbad + 1) {
123  sum += sss_sav[j][i];
124  ct++;
125  }
126  }
127  if (ct > 0)
128  sssref[iy][ix] = sum / ct;
129  }
130  }
131  }
132  } /* end 1-time grid set up */
133 
134  /* locate LL position within reference grid */
135  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), WOASSSNX + 1), 0);
136  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), WOASSSNY + 1), 0);
137 
138  /* compute longitude and latitude of that grid element */
139  xx = i * dx - 180.0 - dx / 2;
140  yy = j * dy - 90.0 - dy / 2;
141 
142  /* bilinearly interpolate, only using full grid boxes */
143  t = (lon - xx) / dx;
144  u = (lat - yy) / dy;
145 
146  if ((sssref[j][i] > sssbad + 1) && (sssref[j][i + 1] > sssbad + 1) &&
147  (sssref[j + 1][i] > sssbad + 1) && (sssref[j + 1][i + 1] > sssbad + 1)) {
148 
149  sss = (1 - t)*(1 - u) * sssref[j ][i ]
150  + t * (1 - u) * sssref[j ][i + 1]
151  + t * u * sssref[j + 1][i + 1]
152  + (1 - t) * u * sssref[j + 1][i ];
153 
154  } else
155  sss = sssbad;
156 
157  return (sss);
158 }
159 
160 /*
161  * for the HYCOM file usage...
162  */
163 #define HYCOMNX 1440
164 #define HYCOMNY 720
165 
166 float get_hycom_sss(char *sssfile, float lon, float lat, float *sss)
167 /*******************************************************************
168 
169  get_hycom_sss
170 
171  purpose: read in the daiy HYCOM sss and give the interpolated value
172  for the input lat, lon
173 
174  Returns type: int - return status: 0 is good, 1 if any other error
175  will fail with -1 status if an h5 io error found (not associated with
176  the file not being h5 or correct format of h5)
177 
178  Parameters: (in calling order)
179  Type Name I/O Description
180  ---- ---- --- -----------
181  char * sssfile I sss file name
182  float lon I longitude: -180 -> 180
183  float lat I latitude: -90 -> 90
184  float * sss O sss found - will be sssbad
185  if over land or otherwise
186  bad value
187 
188  Modification history:
189  Programmer Date Description of change
190  ---------- ---- ---------------------
191  W. Robinson, SAIC 12 Oct 2012 adapted from get_woasssclim in this file
192  W. Robinson, SAIC 16 Oct 2012 this version 2 will interpolate using
193  0 for any missing data
194 
195  Note that the HYCOM file only contains a dataset of salinity, so
196  the size of 1440, 720 probably should be an ID requirement as is the
197  existance of the dataset name 'salinity'
198 
199  *******************************************************************/ {
200  static int firstCall = 1;
201  static int nx = HYCOMNX;
202  static int ny = HYCOMNY;
203  static float dx = 360.0 / HYCOMNX;
204  static float dy = 180.0 / HYCOMNY;
205 
206  typedef float sssref_t[ HYCOMNX + 2 ];
207  static sssref_t* sssref;
208  static sssref_t* sss_sav;
209  //static float sssref[ HYCOMNY + 2 ][ HYCOMNX + 2 ];
210  //static float sss_sav[ HYCOMNY + 2 ][ HYCOMNX + 2 ];
211 
212  int i, j, iter_ct, iter_max;
213  float t, u, xx, yy;
214 
215  *sss = sssbad;
216  if (firstCall) {
217  sssref = (sssref_t*) malloc((HYCOMNY + 2) * sizeof (sssref_t));
218  sss_sav = (sssref_t*) malloc((HYCOMNY + 2) * sizeof (sssref_t));
219  float ssstmp[ny][nx], sum;
220  h5io_str fid, dsid, grpid;
221  int nobj, *types, ndim, dim_siz[10], sto_len, ix, iy, ct;
222  int lymin, lxmin, lymax, lxmax;
223  char **o_names;
224  H5T_class_t class;
225  hid_t native_typ;
226 
227  firstCall = 0;
228 
229  printf("Loading SSS reference from HYCOM file: %s\n", sssfile);
230  printf("\n");
231 
232  /*
233  * check the file to see that it is correct - just return to caller if
234  * not right format, but exit if something that should be performed can't be
235  */
236  if (h5io_openr(sssfile, 0, &fid) != 0)
237  return 1;
238  /* the group may need setting to work, see if '/' will do */
239  if (h5io_set_grp(&fid, "/", &grpid) != 0) {
240  fprintf(stderr, "-E- %s, %d: Failed to set trivial group\n",
241  __FILE__, __LINE__);
242  exit(-1);
243  }
244  if (h5io_grp_contents(&grpid, &nobj, &o_names, &types) != 0) {
245  fprintf(stderr, "-E- %s, %d: Failed to find contents of sss file\n",
246  __FILE__, __LINE__);
247  exit(-1);
248  }
249  /* first name should be salinity */
250  if (strcmp(o_names[0], "salinity") != 0)
251  return 1;
252 
253  /* set to the salinity dataset and get info */
254  if (h5io_set_ds(&fid, "salinity", &dsid) != 0)
255  return 1;
256  if (h5io_info(&dsid, NULL, &class, &native_typ, &ndim, dim_siz,
257  &sto_len) != 0) {
258  fprintf(stderr, "-E- %s, %d: Failed to find info on salinity dataset\n",
259  __FILE__, __LINE__);
260  exit(-1);
261  }
262 
263  if ((ndim != 2) || (dim_siz[0] != 720) || (dim_siz[1] != 1440))
264  return 1;
265  /*
266  * at this point, we have identified the salinity file, just read
267  * the dataset now
268  */
269  if (h5io_rd_ds(&dsid, (void *) ssstmp) != 0) {
270  fprintf(stderr, "-E- %s, %d: Failed to read salinity dataset\n",
271  __FILE__, __LINE__);
272  exit(-1);
273  }
274  /* close and go */
275  h5io_close(&dsid);
276  h5io_close(&fid);
277 
278  /* add wrapping border to simplify interpolation */
279  /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */
280  for (j = 0; j < ny; j++) {
281  for (i = 0; i < nx; i++) {
282  if ((ssstmp[j][i] > 0) && (ssstmp[j][i] <= 1000.))
283  sssref[j + 1][i + 1] = ssstmp[j][i];
284  else
285  sssref[j + 1][i + 1] = sssbad;
286  }
287  sssref[j + 1][0] = sssref[j + 1][nx];
288  sssref[j + 1][nx + 1] = sssref[j + 1][1];
289  }
290  for (i = 0; i < nx + 2; i++) {
291  sssref[0] [i] = sssref[1][i];
292  sssref[ny + 1][i] = sssref[ny][i];
293  }
294  /*
295  * extend data to have good interpolation points at edge = coast
296  */
297  iter_ct = 0;
298  iter_max = 4;
299  do {
300  iter_ct++;
301  memcpy(sss_sav, sssref,
302  (HYCOMNY + 2) * (HYCOMNX + 2) * sizeof ( float));
303  for (iy = 0; iy < ny + 2; iy++) {
304  lymin = (iy == 0) ? 0 : iy - 1;
305  lymax = (iy == ny + 1) ? ny + 1 : iy + 1;
306 
307  for (ix = 0; ix < nx + 2; ix++) {
308  if (sssref[iy][ix] < sssbad + 1) {
309  lxmin = (ix == 0) ? 0 : ix - 1;
310  lxmax = (ix == nx + 1) ? nx + 1 : ix + 1;
311 
312  sum = 0.;
313  ct = 0;
314  for (j = lymin; j < lymax + 1; j++)
315  for (i = lxmin; i < lxmax + 1; i++) {
316  if (sss_sav[j][i] > sssbad + 1) {
317  sum += sss_sav[j][i];
318  ct++;
319  }
320  }
321  if (ct > 0)
322  sssref[iy][ix] = sum / ct;
323  }
324  }
325  }
326  } while (iter_ct < iter_max);
327  /* and end set up */
328  }
329 
330  /* locate LL position within reference grid */
331  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), nx + 1), 0);
332  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), ny + 1), 0);
333 
334  /* compute longitude and latitude of that grid element */
335  xx = i * dx - 180.0 - dx / 2;
336  yy = j * dy - 90.0 - dy / 2;
337 
338  /* bilinearly interpolate, replacing missing (land) values with
339  average of valid values in box */
340  t = (lon - xx) / dx;
341  u = (lat - yy) / dy;
342 
343  if ((sssref[j][i] > sssbad + 1) && (sssref[j][i + 1] > sssbad + 1) &&
344  (sssref[j + 1][i] > sssbad + 1) && (sssref[j + 1][i + 1] > sssbad + 1))
345  *sss = (1 - t)*(1 - u) * sssref[j ][i ]
346  + t * (1 - u) * sssref[j ][i + 1]
347  + t * u * sssref[j + 1][i + 1]
348  + (1 - t) * u * sssref[j + 1][i ];
349  else
350  *sss = sssbad;
351 
352  return 0;
353 }
354 
355 /* ----------------------------------------------------------------------------------- */
356 /* get_sssref() - retrieves reference sea surface salinity . */
357 /* ----------------------------------------------------------------------------------- */
358 #define WOACLIM 1
359 
360 float get_sssref(char *sssfile, float lon, float lat, int day) {
361  float sss = sssbad;
362  int32_t sd_id;
363 
364  if (format < 0) {
365 
366  /* Does the file exist? */
367  if (access(sssfile, F_OK) || access(sssfile, R_OK)) {
368  printf("-E- %s: SSS input file '%s' does not exist or cannot open.\n",
369  __FILE__, sssfile);
370  exit(1);
371  }
372 
373  /* What is it? */
374  sd_id = SDstart(sssfile, DFACC_RDONLY);
375  if (sd_id != FAIL) {
376  /* Format is HDF-like */
377  char title[255] = "";
378  if (SDreadattr(sd_id, SDfindattr(sd_id, "Title"), (VOIDP) title) == 0) {
379  if (strstr(title, "WOA Sea Surface Salinity Monthly Climatology") != NULL) {
380  format = WOACLIM;
381  }
382  } else {
383  printf("-E- %s: unable to read SSS input file %s.\n", __FILE__, sssfile);
384  exit(1);
385  }
386  SDend(sd_id);
387  }
388  else {
389  if (get_hycom_sss(sssfile, lon, lat, &sss) == 0)
390  format = 2;
391  }
392  }
393 
394  switch (format) {
395  case WOACLIM:
396  sss = get_woasssclim(sssfile, lon, lat, day);
397  break;
398  case 2:
399  get_hycom_sss(sssfile, lon, lat, &sss);
400  break;
401  default:
402  printf("-E- %s: unknown SSS input file format for %s.\n", __FILE__, sssfile);
403  exit(1);
404  break;
405  }
406 
407  // assume fresh water if no data
408 
409  if (sss < 0.0) {
410  sss = 0.0;
411  }
412 
413  return (sss);
414 }
415 
416 
MAX
#define MAX(A, B)
Definition: swl0_utils.h:25
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
MIN
#define MIN(x, y)
Definition: rice.h:169
t
data_t t[NROOTS+1]
Definition: decode_rs.h:77
get_woasssclim
float get_woasssclim(char *sssfile, float lon, float lat, int day)
Definition: sssref.c:22
j
int j
Definition: decode_rs.h:73
day
int32_t day
Definition: atrem_corl1v2.h:308
status
int status
Definition: l1_czcs_hdf.c:32
h5io_openr
int h5io_openr(char *file, int opt, h5io_str *id)
Definition: h5io.c:4
FAIL
#define FAIL
Definition: ObpgReadGrid.h:18
NULL
#define NULL
Definition: decode_rs.h:63
h5io_info
int h5io_info(h5io_str *id, char *attr_name, H5T_class_t *class, hid_t *native_typ, int *ndim, int *dim_siz, int *sto_len)
Definition: h5io.c:173
h5io_close
int h5io_close(h5io_str *id)
Definition: h5io.c:115
h5io_set_grp
int h5io_set_grp(h5io_str *id, char *path_name, h5io_str *grp_id)
Definition: h5io.c:369
WOASSSNX
#define WOASSSNX
Definition: sssref.c:19
l12_proto.h
h5io_set_ds
int h5io_set_ds(h5io_str *id, char *path_name, h5io_str *ds_id)
Definition: h5io.c:324
run_local.lat
lat
Definition: run_local.py:69
color_dtdb.ct
list ct
Definition: color_dtdb.py:272
run_local.lon
lon
Definition: run_local.py:69
WOACLIM
#define WOACLIM
Definition: sssref.c:358
h5io_grp_contents
int h5io_grp_contents(h5io_str *id, int *n_obj, char ***o_names, int **types)
Definition: h5io.c:1506
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
get_hycom_sss
float get_hycom_sss(char *sssfile, float lon, float lat, float *sss)
Definition: sssref.c:166
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
u
data_t u
Definition: decode_rs.h:74
rank
Extra metadata that will be written to the HDF4 file l2prod rank
Definition: HOWTO_Add_a_product.txt:80
HYCOMNX
#define HYCOMNX
Definition: sssref.c:163
WOASSSNY
#define WOASSSNY
Definition: sssref.c:20
HYCOMNY
#define HYCOMNY
Definition: sssref.c:164
get_sssref
float get_sssref(char *sssfile, float lon, float lat, int day)
Definition: sssref.c:360
h5io_rd_ds
int h5io_rd_ds(h5io_str *ds_id, void *data)
Definition: h5io.c:505
i
int i
Definition: decode_rs.h:71


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster