OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
get_dem_height.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 
3 /*
4  Exactly replicate the functionality of get_dem_height.f
5  Originally written by: Frederick S. Patt, SAIC GSC, June 1, 2000
6  Translated by: Gwyn F. Fireman, SAIC, July 2013 (with help from f2c)
7  */
8 
9 /* DEM structure (from dem_s.fin) */
10 struct {
11  int32 nrows; // Number of rows in grid
12  int32 neq; // Number of tiles in equatorial row
13  int32 ntiles; // Number of tiles in grid
14  int32 nftiles; // Number of filled tiles
15  int32 isize; // Tile array size
16  int32 itsize; // Array size of nominal area
17  int32 ioff; // Offset of nominal tile in array
18  int16 ncols[121]; // Number of tiles in each row
19  int16 nstart[121]; // Start tile # for each row (0-base)
20  int16 irecno[18338]; // Record no. for tile (0-base; -1 if not filled)
21  int16 iflag[18338]; // Tile flag (0=sea, 1=flag, 2=filled)
22  int16 ihmin[18338]; // Minimum terrain height in tile
23  int16 ihmax[18338]; // Maximum terrain height in tile
24 } dem;
25 
26 /* Function prototypes */
27 float bilin(float *, float *, int32 *, int16 *, int32 *);
28 int initialize_dem(char *, int32 *);
29 
30 /* get_dem_height */
31 int get_dem_height(char *demfile,
32  float *xlon, float *xlat,
33  float *senz, float *sena,
34  float *height) {
35  /* Local variables */
36  int32 status = 0;
37  static float tilat, tilon, tslon, tslat, tansnz, sinsna, cossna;
38  static float x0, y0, x, y, dx, dy, dd, dist, distn, ddist;
39  static float coslat, hgt1, hgt2, los1, los2;
40  static int16 npts, iflag;
41  static int16 i, irow, icol, itile;
42  static int16* tile;
43  static int32 tileid;
44 
45  /* Initialized data */
46  static int firstCall = 1;
47  static int32 idims[3] = {1, 220, 220};
48  static int32 istart[3] = {0, 0, 0};
49  static int32 itilem = -999;
50  static float step = .5f;
51  static float tszmin = .01f;
52  static float snzmax = 56.f;
53  static float remm = 6371000.f; // Earth mean radius in meters
54 
55  /* Check if file is open */
56  if (firstCall) {
57  firstCall = 0;
58 
59  tile = (int16*) allocateMemory(48400 * sizeof (int16), "get_dem_height - tile");
60 
61  /* Open DEM file and fill structure */
62  if (initialize_dem(demfile, &tileid) != 0) {
63  return 1;
64  }
65 
66  /* Compute constants used for calculations */
67  /* Tile size in latitude */
68  tslat = 180.f / (dem.nrows - 1);
69 
70  /* Nominal tile grid point spacing in meters */
71  distn = remm * tslat / (dem.itsize * RADEG);
72 
73  /* Set step size in meters */
74  ddist = step * distn;
75  }
76 
77  /* Compute tile # */
78  irow = (*xlat - 0.0001f + 90.f) / tslat + .5f;
79  icol = (*xlon - 0.0001f + 180.f) / 360.f * dem.ncols[irow];
80  itile = dem.nstart[irow] + icol;
81  iflag = dem.iflag[itile];
82 
83  /* Set nominal values (sea level) */
84  *height = 0.f;
85  dist = 0.f;
86 
87  /* If sea level, return input values */
88  if (iflag == 0) {
89  return 0;
90  }
91 
92  /* Compute needed trig functions of angles */
93  tansnz = tan(*senz / RADEG);
94  sinsna = sin(*sena / RADEG);
95  cossna = cos(*sena / RADEG);
96  coslat = cos(*xlat / RADEG);
97 
98  /* If flat tile */
99  if (iflag == 1) {
100  npts = 2;
101  *height = (float) dem.ihmax[itile];
102  dist = *height * tansnz * remm / (remm + *height);
103 
104  /* Else terrain tile */
105  } else {
106 
107  /* Check for tile already in memory */
108  if (itile != itilem) {
109 
110  /* Read tile array from file */
111  istart[0] = dem.irecno[itile];
112  // printf("Tile %d record %d\n",itile,istart[0]);
113  itilem = itile;
114  if (SDreaddata(tileid, istart, NULL, idims, tile) != 0) {
115  printf("Error reading tile %d record %d\n",
116  itile, dem.irecno[itile]);
117  return 1;
118  }
119 
120  /* Compute longitude width of tile */
121  tslon = 360.f / dem.ncols[irow];
122 
123  /* Compute coordinates of lower-left corner of nominal tile */
124  tilat = (irow - .5f) * tslat - 90.f;
125  tilon = icol * tslon - 180.f;
126  }
127 
128  /* Get terrain heights under line-of-sight */
129 
130  /* Determine maximum number of points needed */
131  npts = dem.ihmax[itile] * tansnz / ddist + 2;
132 
133  /* Determine location of input lon/lat on tile grid */
134  /* (recall that points are bin centers) */
135  x0 = (*xlon - tilon) / tslon * dem.itsize + dem.ioff + .5f;
136  y0 = (*xlat - tilat) / tslat * dem.itsize + dem.ioff + .5f;
137 
138  /* Get interpolated height at input lon/lat */
139  *height = bilin(&x0, &y0, &dem.isize, tile, &status);
140  dist = 0.f;
141 
142  /* Check for zero height and low-slope tile */
143  if (*height == 0.f && iflag == 2) {
144  return 0;
145  }
146  /* Compute grid step size */
147  dx = step * sinsna;
148  dy = step * cossna;
149 
150  /* Correct horizontal step size for latitude and tile size */
151  dx = dx * tslat / (coslat * tslon);
152 
153  /* Find intersection point */
154 
155  /* Check for sensor zenith close to zero */
156  if (tansnz < tszmin) {
157  dist = *height * tansnz;
158 
159  /* Else step downward along line-of-sight */
160  } else {
161 
162  /* Initialize search parameters */
163  hgt1 = 0.f;
164  hgt2 = 0.f;
165  i = npts;
166  dist = i * ddist;
167 
168  /* Compute LOS height with correction for Earth curvature effect */
169  los1 = dist * (dist * tansnz / 2.f + remm)
170  / (remm * tansnz - dist);
171 
172  /* Continue until line-of-sight crosses terrain */
173  while (los1 > hgt1) {
174  --i;
175  los2 = los1;
176  hgt2 = hgt1;
177  dist = i * ddist;
178  x = x0 + i * dx;
179  y = y0 + i * dy;
180  hgt1 = bilin(&x, &y, &dem.isize, tile, &status);
181 
182  /* Check for current point outside tile area and print warning message */
183  if (status != 0 && i < 1) {
184  printf("Warning: %d %d outside tile range\n", (int) x, (int) y);
185  }
186  los1 = dist * (dist * tansnz / 2.f + remm)
187  / (remm * tansnz - dist);
188  }
189 
190  /* Interpolate to find intersection point */
191  dd = (hgt1 - los1) / (hgt1 - los1 + los2 - hgt2);
192  dist += dd * ddist;
193  x += dd * dx;
194  y += dd * dy;
195  hgt1 = bilin(&x, &y, &dem.isize, tile, &status);
196 
197  /* Check for final height outside tile area */
198  if (status != 0) {
199 
200  /* If within GAC swath, return with error */
201  if (*senz <= snzmax) {
202  printf("Bilinear interpolation error %d %d %d\n",
203  (int) x, (int) y, (int) dem.isize);
204  return 0;
205  } else {
206  status = 0;
207  }
208 
209  } else {
210  *height = hgt1;
211  }
212  }
213 
214  /* Compute adjustments to lon, lat, solar zenith and azimuth */
215  /* (need to look at corrections for polar regions) */
216  *xlon += dist * sinsna * RADEG / (remm * coslat);
217  *xlat += dist * cossna * RADEG / remm;
218  *senz -= dist * RADEG / remm;
219 
220  /* Check for longitude transition over date line */
221  if (*xlon > 180.f) {
222  *xlon -= 360.f;
223  }
224  if (*xlon < -180.f) {
225  *xlon += 360.f;
226  }
227  }
228 
229  return 0;
230 } /* get_dem_height */
231 
232 
233 #define tile_ref(a_1,a_2) tile[(a_2)*tile_dim1 + a_1]
234 
235 float bilin(float *x, float *y, int32 *isize, int16 *tile, int32 *status) {
236  /* System generated locals */
237  int tile_dim1, tile_offset;
238  float ret_val;
239 
240  /* Local variables */
241  static float dx, dy;
242  static int32 ix, iy;
243 
244  /* Parameter adjustments */
245  tile_dim1 = *isize;
246  tile_offset = 1 + tile_dim1;
247  tile -= (int32) tile_offset;
248 
249  /* Compute indices */
250  ix = *x;
251  iy = *y;
252 
253  /* Check indices against array size */
254  if ((ix > 0) && (ix < *isize) && (iy > 0) && (iy < *isize)) {
255 
256  dx = *x - ix;
257  dy = *y - iy;
258 
259  ret_val = (1 - dx) * (1 - dy) * tile_ref(ix, iy)
260  + dx * (1 - dy) * tile_ref(ix + 1, iy)
261  + (1 - dx) * dy * tile_ref(ix, iy + 1)
262  + dx * dy * tile_ref(ix + 1, iy + 1);
263  status[0] = 0;
264 
265  } else {
266  ret_val = 0.f;
267  status[0] = -1;
268  }
269 
270  return ret_val;
271 } /* bilin */
272 #undef tile_ref
273 
274 int initialize_dem(char *demfile, int32 *tileid) {
275  /* Initialized data */
276  int32 status = 0;
277  static int32 istart = 0;
278  static int32 istr = 1;
279 
280  /* Local variables */
281  static int32 ind;
282  static int32 at_id, sd_id, idims, fileid;
283 
284  status = 0;
285  tileid[0] = -1;
286 
287  /* Open file */
288  fileid = SDstart(demfile, DFACC_RDONLY);
289  if (fileid == -1) {
290  printf("Error opening DEM file: %s\n", demfile);
291  status = -1;
292  return 0;
293  }
294 
295  /* Read file attributes */
296 
297  at_id = SDfindattr(fileid, "Number of rows");
298  if (at_id == -1) {
299  printf("Error getting index for Number of rows\n");
300  return 0;
301  }
302  status = SDreadattr(fileid, at_id, &dem.nrows);
303  if (status == -1) {
304  printf("Error getting value for Number of rows\n");
305  return 0;
306  }
307 
308  at_id = SDfindattr(fileid, "Equator tiles");
309  if (at_id == -1) {
310  printf("Error getting index for Equator tiles\n");
311  return 0;
312  }
313  status = SDreadattr(fileid, at_id, &dem.neq);
314  if (status == -1) {
315  printf("Error getting value for Equator tiles\n");
316  return 0;
317  }
318 
319  at_id = SDfindattr(fileid, "Array size");
320  if (at_id == -1) {
321  printf("Error getting index for Array size\n");
322  return 0;
323  }
324  status = SDreadattr(fileid, at_id, &dem.isize);
325  if (status == -1) {
326  printf("Error getting value for Array size\n");
327  return 0;
328  }
329 
330  at_id = SDfindattr(fileid, "Tile size");
331  if (at_id == -1) {
332  printf("Error getting index for Tile size\n");
333  return 0;
334  }
335  status = SDreadattr(fileid, at_id, &dem.itsize);
336  if (status == -1) {
337  printf("Error getting value for Tile size\n");
338  return 0;
339  }
340 
341  at_id = SDfindattr(fileid, "Array tile offset");
342  if (at_id == -1) {
343  printf("Error getting index for Array tile offset\n");
344  return 0;
345  }
346  status = SDreadattr(fileid, at_id, &dem.ioff);
347  if (status == -1) {
348  printf("Error getting value for Array tile offset\n");
349  return 0;
350  }
351 
352  at_id = SDfindattr(fileid, "Number of tiles");
353  if (at_id == -1) {
354  printf("Error getting index for Number of tiles\n");
355  return 0;
356  }
357  status = SDreadattr(fileid, at_id, &dem.ntiles);
358  if (status == -1) {
359  printf("Error getting value for Number of tiles\n");
360  return 0;
361  }
362 
363  at_id = SDfindattr(fileid, "Number of filled tiles");
364  if (at_id == -1) {
365  printf("Error getting index for Number of filled tiles\n");
366  return 0;
367  }
368  status = SDreadattr(fileid, at_id, &dem.nftiles);
369  if (status == -1) {
370  printf("Error getting value for Number of filled tiles\n");
371  return 0;
372  }
373 
374  /* Read index arrays */
375 
376  idims = dem.nrows;
377 
378  ind = SDnametoindex(fileid, "Row_number_of_tiles");
379  if (ind == -1) {
380  printf("Error getting index for Row_number_of_tiles\n");
381  return 0;
382  }
383  sd_id = SDselect(fileid, ind);
384  if (sd_id == -1) {
385  printf("Error selecting Row_number_of_tiles\n");
386  return 0;
387  }
388  status = SDreaddata(sd_id, &istart, &istr, &idims, dem.ncols);
389  if (status == -1) {
390  printf("Error reading Row_number_of_tiles\n");
391  return 0;
392  }
393  status = SDendaccess(sd_id);
394 
395  ind = SDnametoindex(fileid, "Row_start_tile");
396  if (ind == -1) {
397  printf("Error getting index for Row_start_tile\n");
398  return 0;
399  }
400  sd_id = SDselect(fileid, ind);
401  if (sd_id == -1) {
402  printf("Error selecting Row_start_tile\n");
403  return 0;
404  }
405  status = SDreaddata(sd_id, &istart, &istr, &idims, dem.nstart);
406  if (status == -1) {
407  printf("Error reading Row_start_tile\n");
408  return 0;
409  }
410  status = SDendaccess(sd_id);
411 
412  idims = dem.ntiles;
413 
414  ind = SDnametoindex(fileid, "DEM_tile_record");
415  if (ind == -1) {
416  printf("Error getting index for DEM_tile_record\n");
417  return 0;
418  }
419  sd_id = SDselect(fileid, ind);
420  if (sd_id == -1) {
421  printf("Error selecting DEM_tile_record\n");
422  return 0;
423  }
424  status = SDreaddata(sd_id, &istart, &istr, &idims, dem.irecno);
425  if (status == -1) {
426  printf("Error reading DEM_tile_record\n");
427  return 0;
428  }
429  status = SDendaccess(sd_id);
430 
431  ind = SDnametoindex(fileid, "DEM_tile_flag");
432  if (ind == -1) {
433  printf("Error getting index for DEM_tile_flag\n");
434  return 0;
435  }
436  sd_id = SDselect(fileid, ind);
437  if (sd_id == -1) {
438  printf("Error selecting DEM_tile_flag\n");
439  return 0;
440  }
441  status = SDreaddata(sd_id, &istart, &istr, &idims, dem.iflag);
442  if (status == -1) {
443  printf("Error reading DEM_tile_flag\n");
444  return 0;
445  }
446  status = SDendaccess(sd_id);
447 
448  ind = SDnametoindex(fileid, "Tile_minimum_height");
449  if (ind == -1) {
450  printf("Error getting index for Tile_minimum_height\n");
451  return 0;
452  }
453  sd_id = SDselect(fileid, ind);
454  if (sd_id == -1) {
455  printf("Error selecting Tile_minimum_height\n");
456  return 0;
457  }
458  status = SDreaddata(sd_id, &istart, &istr, &idims, dem.ihmin);
459  if (status == -1) {
460  printf("Error reading Tile_minimum_height\n");
461  return 0;
462  }
463  status = SDendaccess(sd_id);
464 
465  ind = SDnametoindex(fileid, "Tile_maximum_height");
466  if (ind == -1) {
467  printf("Error getting index for Tile_maximum_height\n");
468  return 0;
469  }
470  sd_id = SDselect(fileid, ind);
471  if (sd_id == -1) {
472  printf("Error selecting Tile_maximum_height\n");
473  return 0;
474  }
475  status = SDreaddata(sd_id, &istart, &istr, &idims, dem.ihmax);
476  if (status == -1) {
477  printf("Error reading Tile_maximum_height\n");
478  return 0;
479  }
480  status = SDendaccess(sd_id);
481 
482  /* Get SDS ID for tile array */
483  ind = SDnametoindex(fileid, "DEM_tile_data");
484  if (ind == -1) {
485  printf("Error getting index for DEM_tile_data\n");
486  return 0;
487  }
488  tileid[0] = SDselect(fileid, ind);
489  if (tileid[0] == -1) {
490  printf("Error selecting DEM_tile_flag\n");
491  }
492  return 0;
493 } /* initialize_dem */
494 
495 /**************************************************************************/
496 
497 int interp_dem_height(char *demfile,
498  float *xlon, float *xlat,
499  float *height) {
500  /* Local variables */
501  int32 status = 0;
502  static float tilat, tilon, tslon, tslat;
503  static float x0, y0;
504  static int16 iflag;
505  static int16 irow, icol, itile;
506  static int16* tile;
507  static int32 tileid;
508 
509  /* Initialized data */
510  static int firstCall = 1;
511  static int32 idims[3] = {1, 220, 220};
512  static int32 istart[3] = {0, 0, 0};
513  static int32 itilem = -999;
514 
515  /* Check if file is open */
516  if (firstCall) {
517  firstCall = 0;
518 
519  tile = (int16*) allocateMemory(48400 * sizeof (int16), "interp_dem_height - tile");
520 
521  /* Open DEM file and fill structure */
522  if (initialize_dem(demfile, &tileid) != 0) {
523  return 1;
524  }
525 
526  /* Compute constants used for calculations */
527  /* Tile size in latitude */
528  tslat = 180.f / (dem.nrows - 1);
529  }
530 
531  /* Compute tile # */
532  irow = (*xlat - 0.0001f + 90.f) / tslat + .5f;
533  icol = (*xlon - 0.0001f + 180.f) / 360.f * dem.ncols[irow];
534  itile = dem.nstart[irow] + icol;
535  iflag = dem.iflag[itile];
536 
537  /* Set nominal values (sea level) */
538  *height = 0.f;
539 
540  /* If sea level, return input values */
541  if (iflag == 0) {
542  return 0;
543  }
544 
545  /* If flat tile */
546  if (iflag == 1) {
547  *height = (float) dem.ihmax[itile];
548 
549  /* Else terrain tile */
550  } else {
551 
552  /* Check for tile already in memory */
553  if (itile != itilem) {
554 
555  /* Read tile array from file */
556  istart[0] = dem.irecno[itile];
557  // printf("Tile %d record %d\n",itile,istart[0]);
558  itilem = itile;
559  if (SDreaddata(tileid, istart, NULL, idims, tile) != 0) {
560  printf("Error reading tile %d record %d\n",
561  itile, dem.irecno[itile]);
562  return 1;
563  }
564 
565  /* Compute longitude width of tile */
566  tslon = 360.f / dem.ncols[irow];
567 
568  /* Compute coordinates of lower-left corner of nominal tile */
569  tilat = (irow - .5f) * tslat - 90.f;
570  tilon = icol * tslon - 180.f;
571  }
572 
573  /* Determine location of input lon/lat on tile grid */
574  /* (recall that points are bin centers) */
575  x0 = (*xlon - tilon) / tslon * dem.itsize + dem.ioff + .5f;
576  y0 = (*xlat - tilat) / tslat * dem.itsize + dem.ioff + .5f;
577 
578  /* Get interpolated height at input lon/lat */
579  *height = bilin(&x0, &y0, &dem.isize, tile, &status);
580  }
581 
582  return 0;
583 }
integer, parameter int16
Definition: cubeio.f90:3
int16 iflag[18338]
int status
Definition: l1_czcs_hdf.c:32
int32 isize
float bilin(float *, float *, int32 *, int16 *, int32 *)
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
int16 ihmin[18338]
#define NULL
Definition: decode_rs.h:63
int32 itsize
int get_dem_height(char *demfile, float *xlon, float *xlat, float *senz, float *sena, float *height)
struct @126 dem
character(len=1000) if
Definition: names.f90:13
int32 ntiles
int32 neq
double precision function f(R1)
Definition: tmd.lp.f:1454
int16 ihmax[18338]
int32 ioff
int32 nrows
int32 nftiles
int initialize_dem(char *, int32 *)
#define RADEG
Definition: czcs_ctl_pt.c:5
#define tile_ref(a_1, a_2)
float xlon[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:93
int16 ncols[121]
int interp_dem_height(char *demfile, float *xlon, float *xlat, float *height)
int16 irecno[18338]
int i
Definition: decode_rs.h:71
int16 nstart[121]