Friday, February 26, 2016

Plot earthquake heatmap on Basemap and Google Map

In the previous blog, we talked about how to query USGS directly online, and plot the earthquakes on the map. In this blog, we will focus on how to use heatmap to represent the earthquakes. We will first try plot the heatmap on Basemap, which is not beautiful, then we will try to generate a nicer figure by plotting on google map.You can find all the script at Qingkai's Github.

Query USGS earthquake catalog

This part is the same as the previous blog.
In [1]:
from bs4 import BeautifulSoup
import urllib
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from functions import *
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [2]:
# let's get the earthquake larger than M5 globally from 2010-01-01 to 2016-01-01. 
url = build_query(outFormat = 'text', starttime = '2010-01-01', endtime = '2016-01-01', minmagnitude = 5.0)
print url

# get the earthquake data from USGS and parse them into a numpy array
r = urllib.urlopen(url).read()
soup = BeautifulSoup(r, "lxml")
events_mat = parse_result(soup.text)
http://earthquake.usgs.gov/fdsnws/event/1/query?format=text&starttime=2010-01-01&endtime=2016-01-01&minmagnitude=5.0
In [3]:
# extract lat and lon for Basemap
lats = [float(item[2]) for item in events_mat]
lons = [float(item[3]) for item in events_mat]

Plot heatmap on Basemap

In [5]:
from matplotlib.colors import LinearSegmentedColormap

m = Basemap(projection='merc',llcrnrlat=-82,urcrnrlat=82,\
            llcrnrlon=-180,urcrnrlon=180, resolution = 'l')
m.drawcoastlines()
m.drawcountries()

# compute appropriate bins to aggregate data
# nx is number of bins in x-axis, i.e. longitude
# ny is number of bins in y-axis, i.e. latitude
nx = 36 # 10 degree for longitude bin
ny = 18 # 10 degree for latitude bin

# form the bins
lon_bins = np.linspace(-180, 180, nx)
lat_bins = np.linspace(-90, 90, ny)
    
# aggregate the number of earthquakes in each bin, we will only use the density
density, lat_edges, lon_edges = np.histogram2d(lats, lons, [lat_bins, lon_bins])

# get the mesh for the lat and lon
lon_bins_2d, lat_bins_2d = np.meshgrid(lon_bins, lat_bins)

# convert the bin mesh to map coordinates:
xs, ys = m(lon_bins_2d, lat_bins_2d) # will be plotted using pcolormesh

# define custom colormap, white -> red, #E6072A = RGB(0.9,0.03,0.16)
cdict = {'red':  ( (0.0,  1.0,  1.0),
                   (1.0,  0.9,  1.0) ),
         'green':( (0.0,  1.0,  1.0),
                   (1.0,  0.03, 0.0) ),
         'blue': ( (0.0,  1.0,  1.0),
                   (1.0,  0.16, 0.0) ) }
custom_map = LinearSegmentedColormap('custom_map', cdict)
plt.register_cmap(cmap=custom_map)

# Here adding one row and column at the end of the matrix, so that 
# density has same dimension as xs, ys, otherwise, using shading='gouraud'
# will raise error
density = np.hstack((density,np.zeros((density.shape[0],1))))
density = np.vstack((density,np.zeros((density.shape[1]))))

# Plot heatmap with the custom color map
plt.pcolormesh(xs, ys, density, cmap="custom_map", shading='gouraud')

# Add color bar and 
cbar = plt.colorbar(orientation='horizontal', shrink=0.625, aspect=20, fraction=0.2,pad=0.02)
cbar.set_label('Number of earthquakes',size=18)

# Plot blue scatter plot of epicenters above the heatmap:    
x,y = m(lons, lats)
m.plot(x, y, 'o', markersize=5,zorder=6, markerfacecolor='#424FA4',markeredgecolor="none", alpha=0.1)

# make image bigger:
plt.gcf().set_size_inches(12,12)

plt.show()

Plot heatmap on Google Map

In [ ]:
!pip install gmaps
In [6]:
import gmaps
In [7]:
# heatmap input data is (lat, lon, weight), weight is optional. 
data = [(float(item[2]), float(item[3])) for item in events_mat]
In [11]:
# plot the heatmap
gmaps.heatmap(data)

Another package to plot on google map

Besides gmaps, another package, which is more flexible, is gmplot, you can use it to generate html directly.
In [ ]:
!pip install gmplot
In [9]:
import gmplot

# declare the center of the map, and how much we want the map zoomed in
gmap = gmplot.GoogleMapPlotter(0, 0, 2)
# plot heatmap
gmap.heatmap(lats, lons)
In [10]:
# save it to html
gmap.draw("Earthquake_heatmap.html")

10 comments:

  1. Very helpful. Thank you Qinghai!

    ReplyDelete
    Replies
    1. You are welcome, glad it helps.

      Delete
    2. Hi Qingkai

      Is lats and longs an array of numbers? I am trying to do exactly this but for different experiment. I have a text file with thousands of latitudes and longitudes in each line. I have read them to arrays like this


      x_data = [44.2, 53.3478, -23.5477, 38.4127, 32.8073, -23.5477]
      y_data = [27.3333, -6.2597, -46.6358, 27.1384, -117.1324, -46.6358]

      x_data is lat, y_data is long.
      Also i have a 2 dim array like [[44.2, 27.3333], [53.3478, -6.2597], [-23.5477, -46.6358], [38.4127, 27.1384], [32.8073, -117.1324], [-23.5477, -46.6358]] which are lat, long pairs

      I used your example, but only gmap.plot works right. I tried with list of 700 GPS, gmap.plot is a very busy image, but all other functions return a blank. Map is present but no dots are represented. What am i doing wrong.

      gmap = gmplot.GoogleMapPlotter(0, 0, 2)

      gmap.plot(x_data, y_data, 'red', edge_width=5)
      #Above API works very well

      # plot heatmap, nothing gets plotted
      #gmap.heatmap(x_data, y_data)

      #below two are blank too
      #gmap.scatter(x_data, y_data, 'r', 100, marker=False)
      #gmap.scatter(x_data, y_data, c='r', marker=True)

      # save it to html
      gmap.draw("country_heatmap.html")

      Thank you!

      Delete
    3. I am not sure what's wrong with your code, but when I am testing, I don't have any issues:

      import gmplot
      import numpy as np

      # generate 700 random lats and lons
      x_data = (np.random.random_sample(size = 700) - 0.5) * 180
      y_data = (np.random.random_sample(size = 700) - 0.5) * 360

      gmap = gmplot.GoogleMapPlotter(0, 0, 2)
      # plot heatmap
      gmap.heatmap(x_data, y_data)
      gmap.scatter(x_data, y_data, c='r', marker=True)

      # save it to html
      gmap.draw("country_heatmap.html")

      Delete
    4. Thanks : i am a new ie to python, so i am sure i am not understanding some areas, thanks to your great help

      My array when i print x_data appears as [27.1384, 27.1384, 27.1384] whereas yours appear as [-30.9088966 57.5068924 70.77143766] with no comma separation.

      I create my array like this
      latitude = float(x)
      longitude = float(y)

      x_data.append(latitude)
      y_data.append(longitude)

      looping it. Initialized x_data = [] and y_data = []

      Also when I run your script with 700, i see lot of green dots, and when i zoom far out i see shades of red. I don't see like your picture on the blog. Does it matter browser/machines etc. I tried both IE and chrome in windows 7. So i see two issues, first i need to see something and then right heatmap..

      thank you very much..

      Delete
    5. No worries.

      If you want to generate the same plot in my blog, you can play with the notebook in the following link:
      https://github.com/qingkaikong/blog/tree/master/10_plot_heatmap_on_Basemap

      Also, the code I sent you before is to generate 700 random dots on the map uniformly, therefore, you can not see the heatmap clear, since they are not clustered together.

      Delete
    6. Thanks can you please let me know your python version, i have 3.4 and upon search i found this
      https://www.reddit.com/r/Python/comments/4bemmb/is_there_a_gmplot_for_python_34/

      i don't understand this fully, so if you can guide me little more much appreciated

      Delete
    7. Ok, I see. I am using python 2.7, all the code was written in 2.7, therefore, it maybe slightly different in version 3.4.

      Delete
  2. Hi I am creating a heatmap using folium and is there anyway to make it so the heatmap doesnt overlay ontop of water?

    ReplyDelete
  3. Five weeks ago my boyfriend broke up with me. It all started when i went to summer camp i was trying to contact him but it was not going through. So when I came back from camp I saw him with a young lady kissing in his bed room, I was frustrated and it gave me a sleepless night. I thought he will come back to apologies but he didn't come for almost three week i was really hurt but i thank Dr.Azuka for all he did i met Dr.Azuka during my search at the internet i decided to contact him on his email dr.azukasolutionhome@gmail.com he brought my boyfriend back to me just within 48 hours i am really happy. What’s app contact : +44 7520 636249‬

    ReplyDelete