Extract raster data value at a point

raster-point.py

from __future__ import division
from osgeo import gdal
from geopy.geocoders import Nominatim

def get_value_at_point(rasterfile, pos):
  gdata = gdal.Open(rasterfile)
  gt = gdata.GetGeoTransform()
  data = gdata.ReadAsArray().astype(np.float)
  gdata = None
  x = int((pos[0] - gt[0])/gt[1])
  y = int((pos[1] - gt[3])/gt[5])
  return data[y, x]

city = 'Ithaca, NY'
g = Nominatim().geocode(city, timeout=5)
p = (g.longitude,g.latitude)
print get_value_at_point('path/to/my/raster/file', p)