Get Elevation for Multiple Coordinates in Python – How to

pandaspython

I have seen the answer to Elevation info for multiple lat long coordinates and I would like to be able to do this with a python script instead. I cannot use Googles API because I have 8000 points. I have them as part of a .csv file which I import into pandas as a dataframe.

something like:

elevation_list = []

def elevation_function(df):
    for lat,long in zip(df.lat,df.lon):
          elevation_i = some_call_to_file(lat,long)
          elevation_list.append(elevation_i)

My coordinates span the western half of the US so Im not sure the function could be so simple. I have downloaded all SRTM90M resolution tiffs from this amazing person

Best Answer

Don't know if this is the simplest way, but it saves gathering elevation data. The USGS-National Map has a REST service that you can use to query elevation for lat/lon coords.

Service url: https://nationalmap.gov/epqs/

You can use pythons requests library and format your query string according to the service parameters. You need your input coordinates in NAD83 (lat/lon).

import requests
import urllib
import pandas as pd

# USGS Elevation Point Query Service
url = r'https://nationalmap.gov/epqs/pqs.php?'

# coordinates with known elevation 
lat = [48.633, 48.733, 45.1947, 45.1962]
lon = [-93.9667, -94.6167, -93.3257, -93.2755]

# create data frame
df = pd.DataFrame({
    'lat': lat,
    'lon': lon
})

def elevation_function(df, lat_column, lon_column):
    """Query service using lat, lon. add the elevation values as a new column."""
    elevations = []
    for lat, lon in zip(df[lat_column], df[lon_column]):

        # define rest query params
        params = {
            'output': 'json',
            'x': lon,
            'y': lat,
            'units': 'Meters'
        }

        # format query string and return query value
        result = requests.get((url + urllib.parse.urlencode(params)))
        elevations.append(result.json()['USGS_Elevation_Point_Query_Service']['Elevation_Query']['Elevation'])

    df['elev_meters'] = elevations

elevation_function(df, 'lat', 'lon')
df.head()

Related Solutions

Python Raster Elevation – How to Get Elevation at Lat/Long Using Python

Here's a more programmatic way of using GDAL than @Aragon's answer. I've not tested it, but it is mostly boiler-plate code that has worked for me in the past. It relies on Numpy and GDAL bindings, but that's about it.

import osgeo.gdal as gdal
import osgeo.osr as osr
import numpy as np
from numpy import ma

def maFromGDAL(filename):
    dataset = gdal.Open(filename, gdal.GA_ReadOnly)

    if dataset is None:
        raise Exception()

    # Get the georeferencing metadata.
    # We don't need to know the CRS unless we want to specify coordinates
    # in a different CRS.
    #projection = dataset.GetProjection()
    geotransform = dataset.GetGeoTransform()

    # We need to know the geographic bounds and resolution of our dataset.
    if geotransform is None:
        dataset = None
        raise Exception()

    # Get the first band.
    band = dataset.GetRasterBand(1)
    # We need to nodata value for our MaskedArray later.
    nodata = band.GetNoDataValue()
    # Load the entire dataset into one numpy array.
    image = band.ReadAsArray(0, 0, band.XSize, band.YSize)
    # Close the dataset.
    dataset = None

    # Create a numpy MaskedArray from our regular numpy array.
    # If we want to be really clever, we could subclass MaskedArray to hold
    # our georeference metadata as well.
    # see here: http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
    # For details.
    masked_image = ma.masked_values(image, nodata, copy=False)
    masked_image.fill_value = nodata

    return masked_image, geotransform

def pixelToMap(gt, pos):
    return (gt[0] + pos[0] * gt[1] + pos[1] * gt[2],
            gt[3] + pos[0] * gt[4] + pos[1] * gt[5])

# Reverses the operation of pixelToMap(), according to:
# https://en.wikipedia.org/wiki/World_file because GDAL's Affine GeoTransform
# uses the same values in the same order as an ESRI world file.
# See: http://www.gdal.org/gdal_datamodel.html
def mapToPixel(gt, pos):
    s = gt[0] * gt[4] - gt[3] * gt[1]
    x = (gt[4] * pos[0] - gt[1] * pos[1] + gt[1] * gt[5] - gt[4] * gt[2]) / s
    y = (-gt[3] * pos[0] + gt[0] * pos[1] + gt[3] * gt[2] - gt[0] * gt[5]) / s
    return (x, y)

def valueAtMapPos(image, gt, pos):
    pp = mapToPixel(gt, pos)
    x = int(pp[0])
    y = int(pp[1])

    if x < 0 or y < 0 or x >= image.shape[1] or y >= image.shape[0]:
        raise Exception()

    # Note how we reference the y column first. This is the way numpy arrays
    # work by default. But GDAL assumes x first.
    return image[y, x]

try:
    image, geotransform = maFromGDAL('myimage.tif')
    val = valueAtMapPos(image, geotransform, (434323.0, 2984745.0))
    print val
except:
    print('Something went wrong.')

[GIS] Assign a point to polygon using pandas and shapely

If you're working with spatial data and Pandas you should take a look at GeoPandas.

The example below demonstrates how to perform a spatial join in GeoPandas (which uses Shapely). A GeoDataFrame object is created from a list of cities and their coordinates and is joined to an ESRI Shapefile containing countries.

import pandas
import geopandas
import geopandas.tools
from shapely.geometry import Point

# Create a DataFrame with some cities, including their location
raw_data = [
    ("London", 51.5, -0.1),
    ("Paris", 48.9, 2.4),
    ("San Francisco", 37.8, -122.4),
]
places = pandas.DataFrame(raw_data, columns=["name", "latitude", "longitude"])

# Create the geometry column from the coordinates
# Remember that longitude is east-west (i.e. X) and latitude is north-south (i.e. Y)
places["geometry"] = places.apply(lambda row: Point(row["longitude"], row["latitude"]), axis=1)
del(places["latitude"], places["longitude"])

# Convert to a GeoDataFrame
places = geopandas.GeoDataFrame(places, geometry="geometry")

# Declare the coordinate system for the places GeoDataFrame
# GeoPandas doesn't do any transformations automatically when performing
# the spatial join. The layers are already in the same CRS (WGS84) so no
# transformation is needed.
places.crs = {"init": "epsg:4326"}

# Load the countries polygons
countries = geopandas.GeoDataFrame.from_file("ne_110m_admin_0_countries.shp")
# Drop all columns except the name and polygon geometry
countries = countries[["name", "geometry"]]

# Perform the spatial join
result = geopandas.tools.sjoin(places, countries, how="left")

# Print the results...
print(result.head())

Note that the spatial join feature is still fairly new and hasn't made it into the stable branch yet - you'll need to download and install the development version.

https://github.com/geopandas/geopandas

The result looks like this:

       name_left             geometry  index_right      name_right
0         London    POINT (-0.1 51.5)           57  United Kingdom
1          Paris     POINT (2.4 48.9)           55          France
2  San Francisco  POINT (-122.4 37.8)          168   United States

You can also use GeoPandas to plot the data with matplotlib:

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1)
countries.plot(ax=ax, color="#cccccc")
places.plot(ax=ax, markersize=5, color="#cc0000")
plt.show()

Best Answer

Related Solutions

Python Raster Elevation – How to Get Elevation at Lat/Long Using Python

[GIS] Assign a point to polygon using pandas and shapely

Related Question