Seismic Events for Q1 2017

Plot it here, there and everywhere...

My friend showed me the Phivolcs' data set on seismic activities and asked if I can interpret the data. I reviewed it and figured it would be great to add to my portfolio. At the same time, my Coursera assignment's due was almost up. We were tasked to pull data from the web and provide interpretation and data visualization. Sort of putting everything we learned all together...being an independent data scientist. I decided to use the Phivolcs data for the assignment. Hitting 2 birds with one stone...heh!

I started to work on the 2014 data since it showed all the data in 1 page. I used BeautifulSoup to scrape the data and I almost went insane after a week of working on the HTML file. I did not get the format that I needed.

So... To make things easier, I just used the 2017 data since the HTML tables were well formatted. I extracted a 3-month data sample of seismic activities from Phivolcs and converted the data to dataframes using Pandas. I used Matplotlib to plot the earthquake magnitudes and leveraged the Basemap toolkit to draw the map and plot the magnitudes in their respective longitude and latitude coordinates. Ha!


The complete python code may be found here. If you interested with the source file, I have extracted the raw logs as and saved as .csv here.

Here is a breakdown of my source code

Import the essential libraries needed. I am running Python 3.5 over Windows 8.1. As far as I can remember, I updated matplotlib to get the basemap library. I ran the command pip install matplotlib update under the administrator account to initiate the task. You may refer here for the documentation on basemap.
from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt import pandas as pd from urllib.request import urlopen

The Q1 data was extracted from Phivolcs' site. The monthly data were encoded on separate pages and were extracted separately. The read_HTML function would have an output of list of dataframes. After inspection, the needed dataframes were on the 3rd index and were extracted respectively. The 3 dataframes were then concatenated and named as df. The dataframe columns were renamed for proper labeling. Saving the file is only optional. I used to_csv to export the file.
#Q1 2017 data only jan_url = "http://www.phivolcs.dost.gov.ph/html/update_SOEPD/EQLatest-Monthly/2017_January.htm" jan_html = urlopen(jan_url).read() jan_table = pd.read_html(jan_html) jan_df=jan_table[3].dropna(how='any') feb_url ="http://www.phivolcs.dost.gov.ph/html/update_SOEPD/EQLatest-Monthly/2017_February.htm" feb_html = urlopen(feb_url).read() feb_table = pd.read_html(feb_html) feb_df=feb_table[3].dropna(how='any') mar_url = "http://www.phivolcs.dost.gov.ph/html/update_SOEPD/EQLatest-Monthly/2017_March.html" mar_html = urlopen(mar_url).read() mar_table = pd.read_html(mar_html) mar_df=mar_table[3].dropna(how='any') df = pd.concat([jan_df,feb_df,mar_df]) df.rename(columns={0:'Date - Time',1:'Latitude',2:'Longitude',3:'Depth (km)',4:'Magnitude',5:'Location' }, inplace=True) #save file df.to_csv('earthquake.csv', sep=',')

Finally! Time to draw! Basemap is pretty easy to follow. (Not!) I used the Transverse Mercator Projection to have a realistic map.I also used the colorbar to aid in the data representation.
#draw map using basemap fig = plt.figure(figsize=(10,10)) plt.title('Philippine Earthquakes \n Q1 of 2017 ') map = Basemap(llcrnrlon=115,llcrnrlat=0,urcrnrlon=130,urcrnrlat=25,projection='tmerc', lat_0 = 12.8797, lon_0 = 121.7740) map.shadedrelief() #plot quake points x, y = map(df.Longitude.tolist(),df.Latitude.tolist()) map.scatter(x, y, c=df.Magnitude,marker='o') map.colorbar().set_label('Magnitude') plt.show()


This tutorial helped me a lot in the visualization using Basemap. I hope it can help you too.

And....the output!

Plotting Temperatures

I recently completed the 2nd assignment in the Coursera course. I had so much fun and had to post it on my blog! Ha!

We were tasked to analyze an 11-year data set which contains the maximum and minimum temperatures for everyday from 2005 to 2015. 2005 to 2014 data was sorted and 2015 data was set aside. The maximum and minimum temperatures per day were plotted. Feb 29 data was removed in the data to keep it clean. Finally, 2015 data was processed and temperature outliers were highlighted.

The following libraries were used.
import pandas as pd import numpy as np import matplotlib.pyplot as plt import matplotlib.patches as mpatches

The data set came from comes from a subset of The National Centers for Environmental Information (NCEI) Daily Global Historical Climatology Network (GHCN-Daily). The GHCN-Daily is comprised of daily climate records from thousands of land surface stations across the globe. 

• The original dataframe has 165085 rows and 4 columns and was assigned to the variable df.
• The 'Data_Value' column is in tenths of degrees C.

In [8]: df = pd.read_csv('course2data.csv') In [9]: df.head() Out[9]: ID Date Element Data_Value 0 USW00094889 11/12/2014 TMAX 22 1 USC00208972 4/29/2009 TMIN 56 2 USC00200032 5/26/2008 TMAX 278 3 USC00205563 11/11/2005 TMAX 139 4 USC00200230 2/27/2014 TMAX -106

Additional columns were added to df for further processing.

• February 29 from leap years were also removed at this stage. 
• 'temp_in_C' was used to convert 'Data_Value' to Celcius. Spyder was able to process the /10 operation but Jupyter kept on crashing. One of the mentors advised that the division operation consumes alot of memory thus causing Jupyter notebook to crash. I opted to use *.10 instead.

#convert date column to real date format df["Date"] = pd.to_datetime(df["Date"]) #spread out month, day and year date_index=pd.DatetimeIndex(df['Date']) df['month']=date_index.month df['day']=date_index.day df['year']=date_index.year df['dayofyear']=date_index.dayofyear # remove feb 29 from dataframe df=df[~((df.month==2) & (df.day ==29))] #add column and convert Temp to celcius #used *.1 instead of /10. division utilizes most of the memory. df['temp_in_C']=df['Data_Value']*.1

Created more dataframes! More fun! At this stage, the maximum and minimum temperatures per day were identified. For example, all temps for November 24 of 2005 to 2014 were collected and the highest (or lowest) temperature was selected and placed to the new dataframe.

• I felt comfortable using pivot_table but groupby function may be used as well.
• 2015 was excluded in the dataframe. We needed to identify which 2015 temperature exceeded the 2005 to 2014 data.

#create df without 2015 data df_data=df[~(df.year==2015)] #create df for TMAX Elements and get max value df_tmax = df_data[df_data.Element=='TMAX'] tmax = pd.pivot_table(df_tmax,index='dayofyear',values=['temp_in_C'],aggfunc={'temp_in_C':'max'}) #create df for TMIN Elements and get min value df_tmin = df_data[df_data.Element=='TMIN'] tmin=pd.pivot_table(df_tmin,index='dayofyear',values=['temp_in_C'],aggfunc={'temp_in_C':'min'}) #create 2015 df and respective tmax and tmin df's df_2015=df[(df.year==2015)] df_2015_tmax = df_2015[df_2015.Element=='TMAX'] tmax_2015=pd.pivot_table(df_2015_tmax,index='dayofyear',values=['temp_in_C'],aggfunc={'temp_in_C':'max'}) df_2015_tmin = df_2015[df_2015.Element=='TMIN'] tmin_2015=pd.pivot_table(df_2015_tmin,index='dayofyear',values=['temp_in_C'],aggfunc={'temp_in_C':'min'}) #create df to show data values that exceed the 2005 to 2014 data tmin_exceed=tmin_2015[tmin_2015.temp_in_C < tmin.temp_in_C.loc[:len(tmin_2015)]] tmax_exceed=tmax_2015[tmax_2015.temp_in_C > tmax.temp_in_C.loc[:len(tmax_2015)]]

Plotting! I no matplotlib expert (yet) but thanks to stackoverflow, I survived. Yey! I did try to minimize the noise to put more emphasis on the data requirement.
#create figure and axes fig = plt.figure() fig, ax = plt.subplots() plt.plot(tmin.temp_in_C,color="white") plt.plot(tmax.temp_in_C,color="white") plt.plot(tmin_exceed.temp_in_C,".",color="blue") plt.plot(tmax_exceed.temp_in_C,".",color="red") # set y axis with degrees celcius softened label by turning to grey ax.yaxis.set_major_formatter(plt.FuncFormatter('{:.0f}$^\circ$C'.format)) plt.yticks(alpha=0.8) #set the range of xticks soften labels by turning to grey plt.xticks(range(0, 366, 30), alpha=0.8) plt.xlabel('Day of Year', alpha=0.8) #set figure title plt.title('Daily High and Low temperatures from 2005-2014 and record breaks of 2015 data.\n',alpha=0.8) # fill the area between the tmax annd tmin data plt.fill_between(range(len(tmax)),tmax.temp_in_C,tmin.temp_in_C,facecolor='blue',alpha=0.25) # remove the frame of the chart for spine in plt.gca().spines.values(): spine.set_visible(False) # remove the ticks (both axes) plt.tick_params(top='off', bottom='on', left='on', right='off', labelleft='on', labelbottom='on') #set up legends red_patch = mpatches.Patch(color='red', label='2015 high temp') blue_patch = mpatches.Patch(color='blue', label='2015 low temp') shade_patch = mpatches.Patch(color='blue', label='2005 to 2014 temp',alpha=0.25) plt.legend(handles=[shade_patch,red_patch,blue_patch],loc=4,frameon=False)

Finally, here is the output.
plt.show()

The complete code may be found here.

TrulyRichClub.com - Do You Want to Gain Financial Wealth and Spiritual Abundance at the Same Time?