# -*- coding: utf-8 -*-

"""

Created on Tue Jan 7 08:48:51 2020

Query NRCAN and create a shapefile using GlobalMapper with the resulting data.

Updates November 10, 2021

Updated time functions to utilize timedeltas

Updated to allow easier change of event type lists

@author: keith.johnson

"""

#!/usr/bin/env python3

import requests as req

import codecs

import os

import datetime

from datetime import date

print("Starting NRCAN script V1.5")

# Sample NRCAN Query from 11/10/2021

# https://www.earthquakescanada.nrcan.gc.ca/fdsnws/event/1/query?endtime=2021-11-08T23%3A59%3A59&eventtype=earthquake%2Cinduced%20or%20triggered%20event&format=text&latitude=63.13&longitude=-90.34&maxdepth=1000&maxlatitude=xx.xxx&maxlongitude=-xxx.xxx&maxmagnitude=10&maxradius=45.04504504504504&mindepth=-5&minlatitude=xx.xxx&minlongitude=-xxx.xxx&minmagnitude=-5&minradius=0&nodata=404&onlyfelt=0&starttime=2021-10-08T00%3A00%3A00

#=======================================================

#Specify Parameters for the Query

#Specify Max/Min Lat Longs

maxlat="xx.xxx"

maxlong="-xxx.xxx"

minlat="xx.xxx"

minlong="-xxx.xxx"

#Aprox area:

#Specify How many days ago to query

daysBack=-30

#quaketypes=L and I, Quakes and Induced

eventtypelist=["earthquake","induced or triggered event"]

#========================================================

#Build event type query from list

# Sample Event type from 11/10/2021

# eventtype=earthquake%2Cinduced%20or%20triggered%20event

eventtype=""

for x in eventtypelist:

eventtype=eventtype+x+","

#format string for URL use

eventtype=eventtype.strip(",")

eventtype=eventtype.replace(" ","%20")

eventtype=eventtype.replace(",","%2C")

eventtype=eventtype+"&"

### DATE SET-UP

#get today's date for file names & date calculations

today=date.today()

Daysbackdelta=datetime.timedelta(days=daysBack)

#Construct starttime and endtime strings

endtime="endtime="+str(today)+"T23%3A59%3A59&"

starttime=str(today+Daysbackdelta)+"T00%3A00%3A00"

#Send request to NRCAN

print("Request Data from NRCAN")

resp = req.get("https://www.earthquakescanada.nrcan.gc.ca/fdsnws/event/1/query?"+eventtype+endtime+"&limit=2000&maxdepth=1000&maxlatitude="+maxlat+"&maxlongitude="+maxlong+"&maxmagnitude=10&mindepth=-5&minlatitude="+minlat+"&minlongitude="+minlong+"&minmagnitude=-5&onlyfelt=0&"+starttime+"&format=text")

#Open file to store NRCAN Response as UTF-8

quakes=codecs.open("query.txt","w","utf-8")

#write UTF-8 to files

quakes.write(resp.text)

#Close file

quakes.close()

#Set outputnames for files

outputname="quakes_"+str(today)+".csv"

shapeoutputname="quakes_"+str(today)+".shp"

print("Converting NRCAN data to CSV")

#Open UTF-8 file and replace | with , delimiters.

with open("query.txt", 'r') as f:

with open(outputname, 'w') as t:

for lines in f:

new_line = lines.replace("|",",")

t.write(new_line)

print("Create Shapefile with Global Mapper")

#Generate GSM file to create shapefile in globalmapper

gsm=open("CSVtoShape.gms","w")

gsm.write("GLOBAL_MAPPER_SCRIPT VERSION=1.00\n")

gsm.write('IMPORT_ASCII FILENAME="'+outputname+'" TYPE=POINT_ONLY COORD_DELIM=COMMA COORD_FORMAT=DECIMAL COORD_ORDER=Y_FIRST SKIP_COLUMNS=2 COL_HEADERS=YES INC_COORD_LINE_ATTRS=YES PROJ=EPSG:4326\n')

gsm.write('EXPORT_VECTOR FILENAME="'+shapeoutputname+'" TYPE=SHAPEFILE SHAPE_TYPE=POINTS\n')

gsm.close()

#Run globalmapper shapefile

os.system('cmd /c C:\Progra~1\GlobalMapper21.0_64bit\global_mapper.exe CSVtoShape.gms')

print("Completed. Written: "+shapeoutputname)