# Data from https://www.ngdc.noaa.gov/mgg/topo/globeget.html
# or google: noaa import GLOBE data into your application
# get the zipfile and the header file using instructions for "Importing into ESRI's ArcGIS 9.x"
# https://www.ngdc.noaa.gov/mgg/topo/report/s11/s11Gxii.html
# headers are here:
# https://www.ngdc.noaa.gov/mgg/topo/elev/esri/arcgis/

import sys
import gdal
import rasterio
FILES = ["e10g.bil", "f10g.bil", "b10g.bil"]

class Elevations():
	def __init__(self):
		pass

	def loadFiles(self):
		self.datasets = {}
		for file in FILES:
			print("loading %s" % file)
			self.datasets[file] = rasterio.open(file)

	def getElevation(self, lon, lat):
		for file in FILES:
			py, px = self.datasets[file].index(lon, lat)
			window = rasterio.windows.Window(px - 1//2, py - 1//2, 1, 1)
			clip = self.datasets[file].read(window=window)
			try:
				return clip[0][0][0]
			except IndexError:
				continue
			return -9999

	def test(self, num):
		v = []
		if num == 1:
			# NS test shows Caribbean, Cayman, Caribbean, Cuba, Gulf of Mexico, US
			lon = -81.25
			for i in range(150, 400):
				lat = i / 10
				e = self.getElevation(lon, lat)
				v.append(e)
			return v
		if num == 2:
			# EW test shows Pacific, Baja, Gulf of Baja, Mexico, Gulf of Mexico, FL, Atlantic
			lat = 28
			for i in range(-1250, -750):
				lon = i / 10
				e = self.getElevation(lon, lat)
				v.append(e)
			return v
		return None

def create():
	ele = Elevations()
	ele.loadFiles()
	return ele

def main():
	ele = create()
	print(ele.test(1))
	print(ele.test(2))

if __name__ == "__main__":
	main()