Search code examples
pythongeocodinggeopandasshapefile

Cant show heighmap on graph

It is not possible to display the graph of the height map. I cant found any problem in code, but graph is empty. Someone can maybe explain how to create shp files for projects like that one. I dont understand: I need to add some fields in shp? Some problem with code?

import matplotlib.pyplot as plt
import geopandas as gpd
import requests
from shapely.geometry import Point
from os.path import join as pth_join
import json
import time

#Configure API request link
request_str = "https://maps.googleapis.com/maps/api/elevation/json?locations=39.7177537,-105.1961914&key=(API here)"

#Define data inputs and outputs
data_root = r"C:\Users\K11El\Desktop\gg"
input_file = "loca.shp"
output_file = input_file[:-4]+"_elev.shp"
#Read input file into a GeoDataFrame
gdf = gpd.read_file(pth_join(data_root,input_file))
gdf_wgs = gdf.to_crs("EPSG:4326")
#Extract x,y from geometry and convert to strings
build_str = gdf_wgs['geometry'].apply(lambda v: ("%s,%s|" % (str(v.y), str(v.x))))
#Limit each request to 512
geom = []
iterations = int(len(build_str)/512) + 1
for x in range(iterations):
    iter_str = build_str[x*512:(x+1)*512]
    #Build concatenated string for API request
    locs = "".join([x for x in iter_str]).rstrip('|')

    #Make request, parse results
    r = requests.get(request_str)
    parsed = json.loads(r.text)

    #Extract elevations from result
    geom = geom + [Point(i['location']['lng'],i['location']
                         ['lat'],i['elevation']) for i in parsed['results']]

    #Slow down API requests so we aren't blocked
    time.sleep(1)
# Save to new file with x,y,z and reproject to CRS of input file
gdf_elev = gpd.GeoDataFrame({'geometry': geom}, crs="EPSG:4326")
gdf_elev.to_crs(gdf.crs).to_file(pth_join(data_root, output_file))
#Check results with a simple plot
gdf_elev['elevation'] = gdf_elev['geometry'].apply(lambda v: (v.z))
gdf_elev.plot("elevation",legend=True,legend_kwds={'label': "Elevation (Meters)"})
plt.show()
Solution
    • using osmnx it's straight forward to use it's implementation of getting elevation from google
    # add x & y columns for osmnx
cities = cities.join(
    cities["geometry"].apply(lambda p: {"x": p.x, "y": p.y}).apply(pd.Series)
)

# create graph
G = ox.utils_graph.graph_from_gdfs(cities, gpd.GeoDataFrame(), graph_attrs=None)
# add elevations using google
G = ox.add_node_elevations_google(G, key)

# now have elevation for all points
gdf = ox.utils_graph.graph_to_gdfs(G, nodes=True, edges=False).set_crs(cities.crs)
    • sample data including elevation
    osmid name_left iso_a3 x y elevation geometry
    0 Vatican City ITA 12.4534 41.9033 24.887 POINT (12.453386544971766 41.903282179960115)
    1 San Marino ITA 12.4418 43.9361 377.216 POINT (12.441770157800141 43.936095834768004)
    2 Rome ITA 12.4813 41.8979 18.527 POINT (12.481312562873995 41.89790148509894)
    3 Vaduz AUT 9.51667 47.1337 513.964 POINT (9.516669472907267 47.13372377429357)
    4 Vienna AUT 16.3647 48.202 185.284 POINT (16.364693096743736 48.20196113681686)
    5 Luxembourg LUX 6.13 49.6117 303.25 POINT (6.130002806227083 49.611660379121076)
    6 Monaco -99 7.40691 43.7396 293.219 POINT (7.406913173465057 43.73964568785249)
    7 Andorra -99 1.51649 42.5 1378.3 POINT (1.5164859605055199 42.5000014435459)
    8 Paris -99 2.33139 48.8686 35.89 POINT (2.33138946713035 48.86863878981461)
    9 Ljubljana SVN 14.515 46.0553 290.773 POINT (14.51496903347413 46.0552883087945)

    full code

    import osmnx as ox
import geopandas as gpd
import pandas as pd
import contextily as ctx

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
cities = gpd.read_file(gpd.datasets.get_path("naturalearth_cities")).sjoin(
    world.loc[world["continent"].eq("Europe")]
)

cities = cities.loc[:, ["name_left", "iso_a3", "geometry"]].reset_index(drop=True)
# add x & y columns for osmnx
cities = cities.join(
    cities["geometry"].apply(lambda p: {"x": p.x, "y": p.y}).apply(pd.Series)
)

# create graph
G = ox.utils_graph.graph_from_gdfs(cities, gpd.GeoDataFrame(), graph_attrs=None)
# add elevations using google
G = ox.add_node_elevations_google(G, key)

# now have elevation for all points
gdf = ox.utils_graph.graph_to_gdfs(G, nodes=True, edges=False).set_crs(cities.crs)


ax = gdf.plot(
    column="elevation",
    figsize=[12, 4],
    markersize=100,
    vmax=gdf["elevation"].quantile(0.9),
    legend=True,
)
gdf.apply(
    lambda x: ax.annotate(
        text=x["name_left"],
        xy=x.geometry.coords[0],
        ha="center",
        xytext=(0, 10),
        textcoords="offset points",
    ),
    axis=1,
)

ctx.add_basemap(ax, crs=gdf.crs, source=ctx.providers.Stamen.Terrain)

    output

    enter image description here