#! /usr/bin/env python

from numpy import *
import numpy as np
import matplotlib.pyplot as plt

# Here we load the topography data file
topo_dat = loadtxt('topogrd.dat')

# Here we define the longitude/latitude grid
lat=180
lon=360
latitude=np.zeros(lat,dtype='f')
longitude=np.zeros(lon,dtype='f')
for i in range(0,lon,1):
   longitude[i]=-180.+i
for j in range(0,lat,1):
   latitude[j]=90.-j

# Here we calculate the topography field on this grid
topo_play=np.zeros((lat,lon),dtype='f')
ii=0
ij=0
for i in range(0,lat,1):
   for j in range(0,lon,1):
      topo_play[i,j]=topo_dat[ii,ij]
      ij=ij+1
      if(ij==10):
         ij=0
         ii=ii+1


# Here we make the histogram
plt.figure(1)
plt.hist(topo_play.flatten(),bins=30)
plt.ylabel('Number of bins')
plt.xlabel('elevation (km)')

# Here we plot the map of the topography
plt.figure(2)
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.contourf(longitude,latitude,topo_play)
plt.colorbar()
plt.show()