3D ski maps. WIP.

hgt.py 1.8KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455
  1. """
  2. This module contains HGT-specific functions
  3. https://gis.stackexchange.com/questions/43743/how-to-extract-elevation-from-hgt-file
  4. These functions are useful for the SRTM dataset which is provided as HGT
  5. files.
  6. """
  7. import struct
  8. def lat_lng_to_row_col(lat, lng):
  9. """
  10. Converts a (latitude, longitude) point to a (row, col) point in the HGT
  11. dataset.
  12. """
  13. # The first row in the file is very likely the northernmost one and there
  14. # are 1200 rows. 3 arc-seconds sampling
  15. return 1201 - round((lng % 1) * 3600 / 3), round((lat % 1) * 3600 / 3)
  16. def find_relevant_points(lat_min, lng_min, lat_max, lng_max, in_file,
  17. verbose=True):
  18. """
  19. Finds the points in the target geographic region which are in in_file.
  20. lat_min, lng_min, lat_max, y_max specify the target geographic
  21. region.
  22. Returns a list of relevant points.
  23. TODO: This function may not be fully working.
  24. """
  25. # Returned list of points
  26. points = []
  27. # Find the corresponding spot in the dataset
  28. row_min, col_min = lat_lng_to_row_col(lat_min, lng_min)
  29. row_max, col_max = lat_lng_to_row_col(lat_max, lng_max)
  30. if col_min > col_max:
  31. col_min, col_max = col_max, col_min
  32. if row_min > row_max:
  33. row_min, row_max = row_max, row_min
  34. with open(in_file, 'rb') as fh:
  35. for i in range(row_max, row_min - 1, -1):
  36. for j in range(col_min, col_max + 1):
  37. # Find the right spot in the file
  38. fh.seek(((i - 1) * 1201 + (j - 1)) * 2)
  39. buf = fh.read(2) # read two bytes and convert them
  40. # ">h" is a signed two byte integer
  41. val = struct.unpack('>h', buf)
  42. # Append (X, Y, Z) to the relevant points list
  43. points.append(j, i, val[0])
  44. return points