[andy/viking.git] / src / dem.c

#include <stdio.h>
#include <string.h>
#include <glib.h>
#include <math.h>
#include <stdlib.h>

#include "dem.h"
#include "file.h"

#define DEM_BLOCK_SIZE 1024
#define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))

static gboolean get_double_and_continue ( gchar **buffer, gdouble *tmp, gboolean warn )
{
  gchar *endptr;
  *tmp = g_strtod(*buffer, &endptr);
  if ( endptr == NULL|| endptr == *buffer ) {
    if ( warn )
      g_warning("Invalid DEM");
    return FALSE;
  }
  *buffer=endptr;
  return TRUE;
}


static gboolean get_int_and_continue ( gchar **buffer, gint *tmp, gboolean warn )
{
  gchar *endptr;
  *tmp = strtol(*buffer, &endptr, 10);
  if ( endptr == NULL|| endptr == *buffer ) {
    if ( warn )
      g_warning("Invalid DEM");
    return FALSE;
  }
  *buffer=endptr;
  return TRUE;
}

static gboolean dem_parse_header ( gchar *buffer, VikDEM *dem )
{
  gdouble val;
  gint int_val;
  guint i;
  gchar *tmp = buffer;

  /* incomplete header */
  if ( strlen(buffer) != 1024 )
    return FALSE;

  /* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
  while (*tmp) {
    if ( *tmp=='D')
      *tmp='E';
    tmp++;
  }

  /* skip name */
  buffer += 149;

  /* "DEM level code, pattern code, palaimetric reference system code" -- skip */
  get_int_and_continue(&buffer, &int_val, TRUE);
  get_int_and_continue(&buffer, &int_val, TRUE);
  get_int_and_continue(&buffer, &int_val, TRUE);

  /* zone */
  get_int_and_continue(&buffer, &int_val, TRUE);
  dem->utm_zone = int_val;
  /* TODO -- southern or northern hemisphere?! */
  dem->utm_letter = 'N';

  /* skip numbers 5-19  */
  for ( i = 0; i < 15; i++ ) {
    if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
      g_warning ("Invalid DEM header");
      return FALSE;
    }
  }

  /* number 20 -- horizontal unit code (utm/ll) */
  get_double_and_continue(&buffer, &val, TRUE);
  dem->horiz_units = val;
  get_double_and_continue(&buffer, &val, TRUE);
  /* dem->orig_vert_units = val; now done below */

  /* TODO: do this for real. these are only for 1:24k and 1:250k USGS */
  if ( dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS ) {
    dem->east_scale = 10.0; /* meters */
    dem->north_scale = 10.0;
    dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
  } else {
    dem->east_scale = 3.0; /* arcseconds */
    dem->north_scale = 3.0;
    dem->orig_vert_units = VIK_DEM_VERT_METERS;
  }

  /* skip next */
  get_double_and_continue(&buffer, &val, TRUE);

  /* now we get the four corner points. record the min and max. */
  get_double_and_continue(&buffer, &val, TRUE);
  dem->min_east = dem->max_east = val;
  get_double_and_continue(&buffer, &val, TRUE);
  dem->min_north = dem->max_north = val;

  for ( i = 0; i < 3; i++ ) {
    get_double_and_continue(&buffer, &val, TRUE);
    if ( val < dem->min_east ) dem->min_east = val;
    if ( val > dem->max_east ) dem->max_east = val;
    get_double_and_continue(&buffer, &val, TRUE);
    if ( val < dem->min_north ) dem->min_north = val;
    if ( val > dem->max_north ) dem->max_north = val;
  }

  return TRUE;
}

static void dem_parse_block_as_cont ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
{
  gint tmp;
  while ( *cur_row < GET_COLUMN(dem, *cur_column)->n_points ) {
    if ( get_int_and_continue(&buffer, &tmp,FALSE) ) {
      if ( dem->orig_vert_units == VIK_DEM_VERT_DECIMETERS )
        GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) (tmp / 10);
      else
        GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) tmp;
    } else
      return;
    (*cur_row)++;
  }
  *cur_row = -1; /* expecting new column */
}

static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
{
  guint n_rows;
  gint i;
  gdouble east_west, south;
  gdouble tmp;

  /* 1 x n_rows 1 east_west south x x x DATA */

  if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
    g_warning("Incorrect DEM Class B record: expected 1");
    return;
  }

  /* don't need this */
  if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;

  /* n_rows */
  if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
    return;
  n_rows = (guint) tmp;

  if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
    g_warning("Incorrect DEM Class B record: expected 1");
    return;
  }
  
  if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
    return;
  if ( !get_double_and_continue(&buffer, &south, TRUE) )
    return;

  /* next three things we don't need */
  if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
  if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
  if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;


  dem->n_columns ++;
  (*cur_column) ++;

  /* empty spaces for things before that were skipped */
  (*cur_row) = (south - dem->min_north) / dem->north_scale;
  if ( south > dem->max_north || (*cur_row) < 0 )
    (*cur_row) = 0;

  n_rows += *cur_row;

  g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
  GET_COLUMN(dem,*cur_column)->east_west = east_west;
  GET_COLUMN(dem,*cur_column)->south = south;
  GET_COLUMN(dem,*cur_column)->n_points = n_rows;
  GET_COLUMN(dem,*cur_column)->points = g_malloc(sizeof(gint16)*n_rows);

  /* no information for things before that */
  for ( i = 0; i < (*cur_row); i++ )
    GET_COLUMN(dem,*cur_column)->points[i] = VIK_DEM_INVALID_ELEVATION;

  /* now just continue */
  dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );


}

static void dem_parse_block ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
{
  /* if haven't read anything or have read all items in a columns and are expecting a new column */
  if ( *cur_column == -1 || *cur_row == -1 ) {
    dem_parse_block_as_header(buffer, dem, cur_column, cur_row);
  } else {
    dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
  }
}

static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename)
{
  gint16 elev;
  gint i, j;
  VikDEM *dem;


  dem = g_malloc(sizeof(VikDEM));

  dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
  dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
  dem->east_scale = dem->north_scale = 3;

  /* TODO */
  dem->min_north = atoi(basename+1) * 3600;
  dem->min_east = atoi(basename+4) * 3600;
  if ( basename[0] == 'S' )
    dem->min_north = - dem->min_north;
  if ( basename[3] == 'W' )
    dem->min_east = - dem->min_east;

  dem->max_north = 3600 + dem->min_north;
  dem->max_east = 3600 + dem->min_east;

  dem->columns = g_ptr_array_new();
  dem->n_columns = 0;

  for ( i = 0; i < 1201; i++ ) {
    dem->n_columns++;
    g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
    GET_COLUMN(dem,i)->east_west = dem->min_east + 3*i;
    GET_COLUMN(dem,i)->south = dem->min_north;
    GET_COLUMN(dem,i)->n_points = 1201;
    GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*1201);
  }

  for ( i = 1200; i >= 0; i-- ) {
    for ( j = 0; j < 1201; j++ ) {
      if ( feof(f) ) {
        g_warning("Incorrect SRTM file: unexpected EOF");
        g_print("%d %d\n", i, j);
        return dem;
      }
      fread(&elev, sizeof(elev), 1, f);
      gchar *x = (gchar *) &elev;
      gchar tmp;
      tmp=x[0];
      x[0]=x[1];
      x[1]=tmp;
      GET_COLUMN(dem,j)->points[i] = elev;
    }

  }

  return dem;
}

#define IS_SRTM_HGT(fn) (strlen((fn))==11 && (fn)[7]=='.' && (fn)[8]=='h' && (fn)[9]=='g' && (fn)[10]=='t' && ((fn)[0]=='N' || (fn)[0]=='S') && ((fn)[3]=='E' || (fn)[3]=='W'))

VikDEM *vik_dem_new_from_file(const gchar *file)
{
  FILE *f;
  VikDEM *rv;
  gchar buffer[DEM_BLOCK_SIZE+1];

  /* use to record state for dem_parse_block */
  gint cur_column = -1;
  gint cur_row = -1;
  const gchar *basename = a_file_basename(file);

      /* FILE IO */
  f = fopen(file, "r");
  if ( !f )
    return NULL;

  if ( strlen(basename)==11 && basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
       (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W'))
    return vik_dem_read_srtm_hgt(f, basename);

      /* Create Structure */
  rv = g_malloc(sizeof(VikDEM));

      /* Header */
  buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
  if ( ! dem_parse_header ( buffer, rv ) ) {
    g_free ( rv );
    return NULL;
  }
  /* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */

  rv->columns = g_ptr_array_new();
  rv->n_columns = 0;

      /* Column -- Data */
  while (! feof(f) ) {
     gchar *tmp;

     /* read block */
     buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';

     /* Fix Fortran-style exponentiation */
     tmp = buffer;
     while (*tmp) {
       if ( *tmp=='D')
         *tmp='E';
       tmp++;
     }

     dem_parse_block(buffer, rv, &cur_column, &cur_row);
  }

     /* TODO - class C records (right now says 'Invalid' and dies) */

  fclose(f);

  /* 24k scale */
  if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
    rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;

  /* FIXME bug in 10m DEM's */
  if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
    rv->min_east -= 100;
    rv->min_north += 200;
  }


  return rv;
}

void vik_dem_free ( VikDEM *dem )
{
  guint i;
  for ( i = 0; i < dem->n_columns; i++)
    g_free ( GET_COLUMN(dem, i)->points );
  g_ptr_array_free ( dem->columns, TRUE );
  g_free ( dem );
}

gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
{
  if ( col < dem->n_columns )
    if ( row < GET_COLUMN(dem, col)->n_points )
      return GET_COLUMN(dem, col)->points[row];
  return VIK_DEM_INVALID_ELEVATION;
}

gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
{
  gint col, row;

  if ( east > dem->max_east || east < dem->min_east ||
      north > dem->max_north || north < dem->min_north )
    return VIK_DEM_INVALID_ELEVATION;

  col = (gint) floor((east - dem->min_east) / dem->east_scale);
  row = (gint) floor((north - dem->min_north) / dem->north_scale);

  return vik_dem_get_xy ( dem, col, row );
}

void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )
{
  *col = (gint) floor((east - dem->min_east) / dem->east_scale);
  *row = (gint) floor((north - dem->min_north) / dem->north_scale);
  if ( *col < 0 ) *col = 0;
  if ( *row < 0 ) *row = 0;
}
Commit	Line	Data
ad0a8c2d EB	1	#include <stdio.h>
	2	#include <string.h>
	3	#include <glib.h>
	4	#include <math.h>
	5	#include <stdlib.h>
	6
	7	#include "dem.h"
	8	#include "file.h"
	9
	10	#define DEM_BLOCK_SIZE 1024
	11	#define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))
	12
	13	static gboolean get_double_and_continue ( gchar *buffer, gdouble tmp, gboolean warn )
	14	{
	15	gchar *endptr;
	16	tmp = g_strtod(buffer, &endptr);
	17	if ( endptr == NULL\|\| endptr == *buffer ) {
	18	if ( warn )
	19	g_warning("Invalid DEM");
	20	return FALSE;
	21	}
	22	*buffer=endptr;
	23	return TRUE;
	24	}
	25
	26
	27	static gboolean get_int_and_continue ( gchar *buffer, gint tmp, gboolean warn )
	28	{
	29	gchar *endptr;
	30	tmp = strtol(buffer, &endptr, 10);
	31	if ( endptr == NULL\|\| endptr == *buffer ) {
	32	if ( warn )
	33	g_warning("Invalid DEM");
	34	return FALSE;
	35	}
	36	*buffer=endptr;
	37	return TRUE;
	38	}
	39
	40	static gboolean dem_parse_header ( gchar buffer, VikDEM dem )
	41	{
	42	gdouble val;
	43	gint int_val;
	44	guint i;
	45	gchar *tmp = buffer;
	46
	47	/* incomplete header */
	48	if ( strlen(buffer) != 1024 )
	49	return FALSE;
	50
	51	/* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
	52	while (*tmp) {
	53	if ( *tmp=='D')
	54	*tmp='E';
	55	tmp++;
	56	}
	57
	58	/* skip name */
	59	buffer += 149;
	60
	61	/* "DEM level code, pattern code, palaimetric reference system code" -- skip */
	62	get_int_and_continue(&buffer, &int_val, TRUE);
	63	get_int_and_continue(&buffer, &int_val, TRUE);
	64	get_int_and_continue(&buffer, &int_val, TRUE);
65
66	/* zone */
67	get_int_and_continue(&buffer, &int_val, TRUE);
68	dem->utm_zone = int_val;
69	/* TODO -- southern or northern hemisphere?! */
70	dem->utm_letter = 'N';
71
72	/* skip numbers 5-19 */
73	for ( i = 0; i < 15; i++ ) {
74	if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
75	g_warning ("Invalid DEM header");
76	return FALSE;
77	}
78	}
79
80	/* number 20 -- horizontal unit code (utm/ll) */
81	get_double_and_continue(&buffer, &val, TRUE);
82	dem->horiz_units = val;
83	get_double_and_continue(&buffer, &val, TRUE);
70147e60	84	/* dem->orig_vert_units = val; now done below */
ad0a8c2d EB	85
	86	/* TODO: do this for real. these are only for 1:24k and 1:250k USGS */
	87	if ( dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS ) {
	88	dem->east_scale = 10.0; /* meters */
	89	dem->north_scale = 10.0;
70147e60	90	dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
ad0a8c2d EB	91	} else {
	92	dem->east_scale = 3.0; /* arcseconds */
	93	dem->north_scale = 3.0;
70147e60	94	dem->orig_vert_units = VIK_DEM_VERT_METERS;
ad0a8c2d EB	95	}
	96
	97	/* skip next */
	98	get_double_and_continue(&buffer, &val, TRUE);
	99
	100	/* now we get the four corner points. record the min and max. */
	101	get_double_and_continue(&buffer, &val, TRUE);
	102	dem->min_east = dem->max_east = val;
	103	get_double_and_continue(&buffer, &val, TRUE);
	104	dem->min_north = dem->max_north = val;
	105
	106	for ( i = 0; i < 3; i++ ) {
	107	get_double_and_continue(&buffer, &val, TRUE);
	108	if ( val < dem->min_east ) dem->min_east = val;
	109	if ( val > dem->max_east ) dem->max_east = val;
	110	get_double_and_continue(&buffer, &val, TRUE);
	111	if ( val < dem->min_north ) dem->min_north = val;
	112	if ( val > dem->max_north ) dem->max_north = val;
	113	}
	114
	115	return TRUE;
	116	}
	117
	118	static void dem_parse_block_as_cont ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	119	{
	120	gint tmp;
	121	while ( cur_row < GET_COLUMN(dem, cur_column)->n_points ) {
70147e60 EB	122	if ( get_int_and_continue(&buffer, &tmp,FALSE) ) {
	123	if ( dem->orig_vert_units == VIK_DEM_VERT_DECIMETERS )
	124	GET_COLUMN(dem, cur_column)->points[cur_row] = (gint16) (tmp / 10);
	125	else
	126	GET_COLUMN(dem, cur_column)->points[cur_row] = (gint16) tmp;
	127	} else
ad0a8c2d EB	128	return;
	129	(*cur_row)++;
	130	}
	131	cur_row = -1; / expecting new column */
	132	}
	133
	134	static void dem_parse_block_as_header ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	135	{
	136	guint n_rows;
	137	gint i;
	138	gdouble east_west, south;
	139	gdouble tmp;
	140
	141	/* 1 x n_rows 1 east_west south x x x DATA */
	142
	143	if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) \|\| tmp != 1 ) {
	144	g_warning("Incorrect DEM Class B record: expected 1");
	145	return;
	146	}
	147
	148	/* don't need this */
	149	if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;
	150
	151	/* n_rows */
	152	if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
	153	return;
	154	n_rows = (guint) tmp;
	155
	156	if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) \|\| tmp != 1 ) {
	157	g_warning("Incorrect DEM Class B record: expected 1");
	158	return;
	159	}
	160
	161	if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
	162	return;
	163	if ( !get_double_and_continue(&buffer, &south, TRUE) )
	164	return;
	165
	166	/* next three things we don't need */
	167	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	168	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	169	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	170
	171
	172	dem->n_columns ++;
	173	(*cur_column) ++;
	174
	175	/* empty spaces for things before that were skipped */
	176	(*cur_row) = (south - dem->min_north) / dem->north_scale;
	177	if ( south > dem->max_north \|\| (*cur_row) < 0 )
	178	(*cur_row) = 0;
	179
	180	n_rows += *cur_row;
	181
	182	g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
	183	GET_COLUMN(dem,*cur_column)->east_west = east_west;
	184	GET_COLUMN(dem,*cur_column)->south = south;
	185	GET_COLUMN(dem,*cur_column)->n_points = n_rows;
	186	GET_COLUMN(dem,cur_column)->points = g_malloc(sizeof(gint16)n_rows);
	187
	188	/* no information for things before that */
	189	for ( i = 0; i < (*cur_row); i++ )
	190	GET_COLUMN(dem,*cur_column)->points[i] = VIK_DEM_INVALID_ELEVATION;
	191
192	/* now just continue */
193	dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );
194
195
196	}
197
198	static void dem_parse_block ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
199	{
200	/* if haven't read anything or have read all items in a columns and are expecting a new column */
201	if ( cur_column == -1 \|\| cur_row == -1 ) {
202	dem_parse_block_as_header(buffer, dem, cur_column, cur_row);
203	} else {
204	dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
205	}
206	}
207
208	static VikDEM vik_dem_read_srtm_hgt(FILE f, const gchar *basename)
209	{
210	gint16 elev;
211	gint i, j;
212	VikDEM *dem;
213
214
215	dem = g_malloc(sizeof(VikDEM));
216
217	dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
70147e60	218	dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
ad0a8c2d EB	219	dem->east_scale = dem->north_scale = 3;
	220
	221	/* TODO */
	222	dem->min_north = atoi(basename+1) * 3600;
	223	dem->min_east = atoi(basename+4) * 3600;
	224	if ( basename[0] == 'S' )
	225	dem->min_north = - dem->min_north;
	226	if ( basename[3] == 'W' )
	227	dem->min_east = - dem->min_east;
	228
	229	dem->max_north = 3600 + dem->min_north;
	230	dem->max_east = 3600 + dem->min_east;
	231
	232	dem->columns = g_ptr_array_new();
	233	dem->n_columns = 0;
	234
	235	for ( i = 0; i < 1201; i++ ) {
	236	dem->n_columns++;
	237	g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
	238	GET_COLUMN(dem,i)->east_west = dem->min_east + 3*i;
	239	GET_COLUMN(dem,i)->south = dem->min_north;
	240	GET_COLUMN(dem,i)->n_points = 1201;
	241	GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*1201);
	242	}
	243
	244	for ( i = 1200; i >= 0; i-- ) {
	245	for ( j = 0; j < 1201; j++ ) {
	246	if ( feof(f) ) {
	247	g_warning("Incorrect SRTM file: unexpected EOF");
	248	g_print("%d %d\n", i, j);
	249	return dem;
	250	}
	251	fread(&elev, sizeof(elev), 1, f);
	252	gchar x = (gchar ) &elev;
	253	gchar tmp;
	254	tmp=x[0];
	255	x[0]=x[1];
	256	x[1]=tmp;
	257	GET_COLUMN(dem,j)->points[i] = elev;
	258	}
	259
	260	}
	261
	262	return dem;
	263	}
	264
	265	#define IS_SRTM_HGT(fn) (strlen((fn))==11 && (fn)[7]=='.' && (fn)[8]=='h' && (fn)[9]=='g' && (fn)[10]=='t' && ((fn)[0]=='N' \|\| (fn)[0]=='S') && ((fn)[3]=='E' \|\| (fn)[3]=='W'))
	266
	267	VikDEM vik_dem_new_from_file(const gchar file)
	268	{
	269	FILE *f;
	270	VikDEM *rv;
	271	gchar buffer[DEM_BLOCK_SIZE+1];
	272
	273	/* use to record state for dem_parse_block */
	274	gint cur_column = -1;
	275	gint cur_row = -1;
	276	const gchar *basename = a_file_basename(file);
	277
	278	/* FILE IO */
	279	f = fopen(file, "r");
	280	if ( !f )
	281	return NULL;
	282
283	if ( strlen(basename)==11 && basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
284	(basename[0] == 'N' \|\| basename[0] == 'S') && (basename[3] == 'E' \|\| basename[3] =='W'))
285	return vik_dem_read_srtm_hgt(f, basename);
286
287	/* Create Structure */
288	rv = g_malloc(sizeof(VikDEM));
289
290	/* Header */
291	buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
292	if ( ! dem_parse_header ( buffer, rv ) ) {
293	g_free ( rv );
294	return NULL;
295	}
296	/* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */
297
298	rv->columns = g_ptr_array_new();
299	rv->n_columns = 0;
300
301	/* Column -- Data */
302	while (! feof(f) ) {
303	gchar *tmp;
304
305	/* read block */
306	buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
307
308	/* Fix Fortran-style exponentiation */
309	tmp = buffer;
310	while (*tmp) {
311	if ( *tmp=='D')
312	*tmp='E';
313	tmp++;
314	}
315
316	dem_parse_block(buffer, rv, &cur_column, &cur_row);
317	}
318
319	/* TODO - class C records (right now says 'Invalid' and dies) */
320
321	fclose(f);
322
323	/* 24k scale */
324	if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
325	rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;
326
327	/* FIXME bug in 10m DEM's */
328	if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
329	rv->min_east -= 100;
330	rv->min_north += 200;
331	}
332
333
334	return rv;
335	}
336
337	void vik_dem_free ( VikDEM *dem )
338	{
339	guint i;
340	for ( i = 0; i < dem->n_columns; i++)
341	g_free ( GET_COLUMN(dem, i)->points );
342	g_ptr_array_free ( dem->columns, TRUE );
343	g_free ( dem );
344	}
345
346	gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
347	{
348	if ( col < dem->n_columns )
349	if ( row < GET_COLUMN(dem, col)->n_points )
350	return GET_COLUMN(dem, col)->points[row];
351	return VIK_DEM_INVALID_ELEVATION;
352	}
353
354	gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
355	{
356	gint col, row;
357
358	if ( east > dem->max_east \|\| east < dem->min_east \|\|
359	north > dem->max_north \|\| north < dem->min_north )
360	return VIK_DEM_INVALID_ELEVATION;
361
362	col = (gint) floor((east - dem->min_east) / dem->east_scale);
363	row = (gint) floor((north - dem->min_north) / dem->north_scale);
364
365	return vik_dem_get_xy ( dem, col, row );
366	}
367
368	void vik_dem_east_north_to_xy ( VikDEM dem, gdouble east, gdouble north, guint col, guint *row )
369	{
370	*col = (gint) floor((east - dem->min_east) / dem->east_scale);
371	*row = (gint) floor((north - dem->min_north) / dem->north_scale);
372	if ( col < 0 ) col = 0;
373	if ( row < 0 ) row = 0;
374	}
375