#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <fcntl.h>
#include <math.h>

/*
8/5/2008

NOTES:
It expects all the SRTM files for a particular resolution to be in a directory;
here we've got /opt/srtm-1arcsecond (with 1117 files) and and /opt/srtm-3arcsecond (14548 files).

Poke around http://www2.jpl.nasa.gov/srtm/ to find the docs for the HGT format

Be careful about sign extension, particularly when comparing the returned elevation
against the -32768 they use to indicate data voids

Compiling on my linux box (the '-std=c99' is because of 'round()'):

  gcc -std=c99 -o readsrtm readsrtm.c -lm

*/


typedef short INT16;

#define BYTES_PER_SAMPLE 2

     /* lat and lon are location whose elevation you want
        dir is where to look for HGT files with the names like S10W010.hgt
        n_samples is number of rows and columns for files in that directory (1201 for 3 arcsecond data and 3601 for 1 arcsecond)
        places the actual data point used in *_lat, *_lon and the elevation in meters in *_elev
        returns 1 on success
        DOESN'T check if the data is "void", you have to do that by comparing result against -32768
    */

int get_elev(double lat, double lon, char *dir, int n_samples, double *_lat, double *_lon, INT16 *_elev) {

     /* ilat0 and ilon0 are the SW corner of the file we want
        The data is packed with the first row of data in the file
        running from the NW to the NE corner, then the next lower latitude, ...
        Adjacent files overlap one row or column, so for instance for the 1 arcsecond files
        the grid spacing is 1/3600 of a degree but there are 3601 samples  */

  int ilat0;
  int ilon0;
  int samples_per_degree;
  int row;
  int col;
  char buf[256];
  unsigned char dbuf[BYTES_PER_SAMPLE + 1];
  int fd;


        /* file is named by it's SW corner */
  ilat0 = floor(lat);
  ilon0 = floor(lon);

  snprintf(buf, sizeof buf, "%s/%c%02d%c%03d.hgt",
				  	  	   dir,
				  	  	   lat >= 0? 'N' : 'S',
				  	  	   abs(ilat0),
				  	  	   lon >= 0? 'E' : 'W',
				  	  	   abs(ilon0));
  if ((fd = open(buf, O_RDONLY)) == -1) {
    perror("get_elev: open failed");
    return 0;
  }


	  /* assumes one row and column of overlap (e.g. grid spacing is 1/3600 degree but there are 3601 samples) */
  samples_per_degree = n_samples - 1;

	  /* find row/column from lat/lon */
	  /* this rounds properly because we know row and col are >= 0 */
  row = round((ilat0 + 1 - lat  ) * samples_per_degree);
  col = round((lon       - ilon0) * samples_per_degree);

  if (lseek(fd, BYTES_PER_SAMPLE * (row * n_samples + col), SEEK_SET) == -1) {
    perror("get_elev: lseek failed");
    close(fd);
    return 0;
  }
  if (read(fd, dbuf, BYTES_PER_SAMPLE) != BYTES_PER_SAMPLE) {
    perror("get_elev: read failed");
    close(fd);
    return 0;
  }

	  /* WARNING: assumes SRTM file is big-endian and we're on a little-endian architecture */
  *_elev = dbuf[1] | (dbuf[0] << 8);


	  /* find lat/lon from row/col (i.e. these are the coordinates of the closest data point) */
  *_lat = ilat0 + 1 - 1. * row / samples_per_degree;
  *_lon = ilon0     + 1. * col / samples_per_degree;

  close(fd);
  return 1;
}



int main(int argc, char *argv[]) {
  double desired_lat;
  double desired_lon;
  double found_lat;
  double found_lon;
  INT16 elev;
  if (argc != 3) {
    fprintf(stderr, "Usage: %s <lat> <lon>\n", argv[0]);
    exit(-1);
  }
  desired_lat = atof(argv[1]);
  desired_lon = atof(argv[2]);

  if (get_elev(desired_lat, desired_lon, "/opt/srtm-1arcsecond", 3601, &found_lat, &found_lon, &elev))
    printf("sought %.5lf %.5lf found %d at %.5lf %.5lf (1\" data)\n", desired_lat, desired_lon, elev, found_lat, found_lon);
  if (get_elev(desired_lat, desired_lon, "/opt/srtm-3arcsecond", 1201, &found_lat, &found_lon, &elev))
    printf("sought %.5lf %.5lf found %d at %.5lf %.5lf (3\" data)\n", desired_lat, desired_lon, elev, found_lat, found_lon);
}