Commit 9d28abc1 authored by Franco Mele's avatar Franco Mele
Browse files

minimalia

parent db21cc74
*.swp
*.history
*.DS_Store
*.ps
# GMT common arguments shelf
-B1::/: :
-Jm0.9
-R4/34/24/49r
-X3.0
-Y1.5
-jm0.9
EOF
EOF
OF
This diff is collapsed.
#x y texti
19.8 48.7 0 1 10 BC Mag
20.5 48.0 0 1 10 BC 6
20.5 47.3 0 1 10 BC 5
20.5 46.8 0 1 10 BC 4
20.5 46.4 0 1 10 BC 3
20.5 46.0 0 1 10 BC 2
17.7 48.7 0 1 10 BC Depth (km)
17.8 48.1 0 1 10 BC ( 0; 5]
17.8 47.65 0 1 10 BC ( 5; 10]
17.9 47.20 0 1 10 BC (10; 30]
17.9 46.75 0 1 10 BC (30; 60]
17.9 46.30 0 1 10 BC (60;100]
17.9 45.85 0 1 10 BC >100
21.8 48.3 0 1 10 BC best
21.8 45.1 0 1 10 BC worst
21.8 48.7 0 1 10 BC Quality
......@@ -21,14 +21,14 @@ GMT_CPT_FILE_MAP1=$3
# Legend file
GMT_LEGEND_CIRCLE=legend_circle_quality
GMT_LEGEND_TEXT=legend_quality
GMT_LEGEND_TEXT=legend_quality4
GMT_TITLE=' '
# Data about original (Italy) map
ITALY_MAP=italy.ps
#ITALY_GRD=italy.grd
#ZONE='6.4565/18.948167/36.670833/47.329167'
ZONE='4/23/34/49'
BAR_POSITION_1='15.75c/15.0c/3.5c/.35c'
ZONE='4/34/24/49r'
BAR_POSITION_1='15.75c/15.0c/5.5c/.35c'
#BAR_POSITION_2='.6c/4.9c/4c/.5c'
# posizione X e Y del foglio in inch.
......@@ -50,7 +50,8 @@ SCAL=0.9
psscale -D$BAR_POSITION_1 -O -C$GMT_CPT_FILE_MAP1 -B1::/:" ": -K >> $ITALY_MAP
# psscale -D$BAR_POSITION_2 -O -C$GMT_CPT_FILE_MAP2 -Ba200f100::/:: -K >> $ITALY_MAP
psxy -R -Jm$SCAL -O $GMT_LEGEND_CIRCLE -S -W0.25p -K >> $ITALY_MAP
pstext -F+f10p $GMT_LEGEND_TEXT -O -R -Jm$SCAL -K >> $ITALY_MAP
# pstext -F+a0+c1+f10 $GMT_LEGEND_TEXT -O -R$ZONE -Jm$SCAL -K >> $ITALY_MAP
pstext $GMT_LEGEND_TEXT -O -R$ZONE -Jm$SCAL -K >> $ITALY_MAP
# Deleting tmp files
#rm $RESAMPLED $PALETTE
......
-100 white 0 253.48/250.43/225.29
0 253.48/250.43/225.29 100 252.48/245.86/195.57
100 252.48/245.86/195.57 200 252/241.57/169.86
200 252/241.57/169.86 300 251.9/237.81/150.76
300 251.9/237.81/150.76 400 250.38/234.76/130.95
400 250.38/234.76/130.95 500 250/231.71/111.14
500 250/231.71/111.14 600 232.67/212.67/93.333
600 232.67/212.67/93.333 700 193.05/173.05/78.095
700 193.05/173.05/78.095 800 153.43/133.43/62.857
800 153.43/133.43/62.857 900 108.1/104.76/50
900 108.1/104.76/50 1000 235.14/253.1/238.24
1000 235.14/253.1/238.24 1100 205.43/249.29/204.86
1100 205.43/249.29/204.86 1200 175.71/245.48/172.1
1200 175.71/245.48/172.1 1300 114.67/213/178.67
1300 114.67/213/178.67 1400 49.143/176.86/190.86
1400 49.143/176.86/190.86 1500 0/136.67/185.71
1500 0/136.67/185.71 1600 0/83.333/128.57
1600 0/83.333/128.57 1700 0/30/71.429
1700 0/30/71.429 1800 0/7.619/38.095
1800 0/7.619/38.095 1900 0/3.8095/19.048
1900 0/3.8095/19.048 2000 black
B white
F black
N white
# cpt file created by: makecpt -I -Crelief -Z -T-100/2000/100
#COLOR_MODEL = RGB
#
-100 255 255 255 0 253 250 225
0 253 250 225 100 252 246 196
100 252 246 196 200 252 242 170
200 252 242 170 300 252 238 151
300 252 238 151 400 250 235 131
400 250 235 131 500 250 232 111
500 250 232 111 600 233 213 93
600 233 213 93 700 193 173 78
700 193 173 78 800 153 133 63
800 153 133 63 900 108 105 50
900 108 105 50 1000 235 253 238
1000 235 253 238 1100 205 249 205
1100 205 249 205 1200 176 245 172
1200 176 245 172 1300 115 213 179
1300 115 213 179 1400 49 177 191
1400 49 177 191 1500 0 137 186
1500 0 137 186 1600 0 83 129
1600 0 83 129 1700 0 30 71
1700 0 30 71 1800 0 8 38
1800 0 8 38 1900 0 4 19
1900 0 4 19 2000 0 0 0
B 255 255 255
F 0 0 0
N 255 255 255
# programmed by Franco Mele on jun 8th, 2021
#
#
#
# programma per l'estrazione di eventi e parametri di qualità da server web (in xml);
# produce il file di input per una procedura che plotta una mappa con le qualità
# delle localizzazioni.
# I criteri per definire le qualità sono
# (vedere Amato - Mele 2008 Performance of the INGV National Seismic Network from 1997 to 2007)
# The first quality rating is based on errors and goodness-of-fit:
# Score
# 1.5 A. RMS < 0.45 sec and ERH <2.0 km and ERZ <4.0 km
# 0.5 B. RMS < 0.90 sec and ERH <5.0 km and ERZ <10.0 km
# -0.5 C. RMS < 1.50 sec and ERH <10.0 km
# -1.5 D. Worse than above
#
# The second quality rating is based on station geometry:
# Score
# 3.0 A. NWR > 6 and MAXGAP < 90 and either DMIN < DEPTH or DMIN <10.0
# 1.0 B. NWR > 6 and MAXGAP < 135 and either DMIN < 2*DEPTH or DMIN <20
# -1.0 C. NWR > 6 and MAXGAP < 180 and DMIN <100
# -3.0 D. Worse than above
# La somma dei due valori corrispondenti alle due qualità può assumere
# 10 valori (-4.5, -3.5, ..-0.5, 0.5, ..., 4.5)
#
# l'uscita del programma è un file di input per GMT contenete 5 colonne:
# longitude , latitude , score_di_qualità , dimensione(M*M/40) , codice-formato(funzione della depth)
#
# 14.0678 43.0577 -2.5 0.24649 n
# 14.3110 38.5290 -2.5 0.18225 n
# 11.1080 46.7760 -3.5 0.24806 a
# 10.2478 44.6922 -1.5 0.15876 s
# 15.2043 38.3988 -1.5 0.20449 a
# 15.7208 39.2178 1.5 0.28730 a
# 15.7802 39.5208 -1.5 0.20880 a
# 14.2312 37.8287 3.5 0.16002 s
#
# i codici GMT sono t, s, n, h, g, a (per i limiti di depth: 0 -> 5, 10, 30, 60, 100, > 100)
# CRITERI ORIGINALI di Hypoinverse2000:
#
# The first quality rating is based on errors and goodness-of-fit:
#
# A. RMS < 0.15 sec and ERH <1.0 km and ERZ <2.0 km
# B. RMS < 0.30 sec and ERH <2.5 km and ERZ <5.0 km
# C. RMS < 0.50 sec and ERH <5.0 km
# D. Worse than above
#
# The second quality rating is based on station geometry:
#
# A. NWR > 6 and MAXGAP < 90 and either DMIN < DEPTH or DMIN <5.0
# B. NWR > 6 and MAXGAP < 135 and either DMIN < 2*DEPTH or DMIN <10
# C. NWR > 6 and MAXGAP < 180 and DMIN <100
# D. Worse than above
import pandas as pd
import requests
import xmltodict
import math
def main():
# http://terremoti.ingv.it/events?last_nd=-1&starttime=2021-05-08&endtime=2021-05-15&minmag=2&maxmag=10&wheretype=pointradius&box_search=Mondo&minlat=-90&maxlat=90&minlon=-180&maxlon=180&municipio=&lat=42&lon=13&maxradiuskm=30&mindepth=-10&maxdepth=1000
typeq = input(" selezione rettangolare (1) o circolare (2) ? ")
if( typeq == "1" ):
minlat = input(" latitudine minima (ES: 34.0) ")
maxlat = input(" latitudine massima (ES: 49.0) ")
minlon = input(" longitudine minima (ES: 4.0) ")
maxlon = input(" longitudine massima (ES: 23.0) ")
elif( typeq =="2" ):
lat = input(" latitudine del centro ")
lon = input(" longitudine del centro ")
radius = input(" raggio del cerchio (km) ")
else:
print(" non hai scelto né 1 né 2 ")
exit()
mindepth = input(" profondita' minima (ES: 0) ")
maxdepth = input(" profondita' massima (ES: 700.0) ")
minmag = input(" magnitudo minima (ES: 2.5) ")
maxmag = input(" magnitudo massima (ES: 10.0) ")
startt = input(" data di inizio nella forma yyyy-mm-dd ")
endt = input(" data di fine nella forma yyyy-mm-dd ")
print(" estraggo i dati del periodo ",startt, " ", endt)
if( typeq =="1" ):
url = "http://webservices.ingv.it/fdsnws/event/1/query?minlat="+minlat+"&maxlat="+maxlat+ \
"&minlon="+minlon+"&maxlon="+maxlon+ \
"&mindepth="+mindepth+"&maxdepth="+maxdepth+ \
"&starttime="+startt+"&endtime="+endt+ \
"&minmag="+minmag+"&maxmag="+maxmag
else:
url = "http://webservices.ingv.it/fdsnws/event/1/query?wheretype=pointradius"+ \
"&lat="+lat+ \
"&lon="+lon+ \
"&maxradiuskm="+radius+ \
"&mindepth="+mindepth+"&maxdepth="+maxdepth+ \
"&starttime="+startt+"&endtime="+endt+ \
"&minmag="+minmag+"&maxmag="+maxmag
url0 = url.replace(" ","")
print(url0)
response = requests.get(url0)
with open('my_text.xml', 'wb') as fil:
fil.write(response.content)
#opening the xml file in read mode
with open("my_text.xml","r") as xml_obj:
#coverting the xml data to Python dictionary
my_dict = xmltodict.parse(xml_obj.read())
#closing the file
xml_obj.close()
#for val in my_dict.values():
# print(type(val))
#for key in my_dict.keys():
# print(key)
#for key in my_dict["q:quakeml"].keys():
# print(key)
#for key in my_dict["q:quakeml"].keys():
# print(key)
# print(" --- START of keys in eventParameter --- ")
# for key in my_dict["q:quakeml"]["eventParameters"].keys():
# print(key)
# print(" --- END of keys in eventParameter --- ")
# print(" --- START of keys in event --- ")
# for val in my_dict["q:quakeml"]["eventParameters"]["event"]:
# print(val)
# print(" --- END of keys in event --- ")
#for val in my_dict["q:quakeml"]["eventParameters"].values():
# print(val, "\n")
filout = startt+"_"+endt+"_quality.txt"
filo = filout.replace("-","")
filo = filo.replace(" ","")
filo = filo.replace(":","")
with open(filo, 'w') as f:
i=0
for event in my_dict["q:quakeml"]["eventParameters"]["event"]:
# print(event)
print(event["description"]["text"])
print(event["origin"]["time"])
if(event["origin"]["creationInfo"].get("version")):
version = event["origin"]["creationInfo"]["version"]
print(" version: ", version)
# ======================================================================
if(version == "1000"):
# ....................=============================================
lat = float(event["origin"]["latitude"]["value"])
print(" lat. ",lat," cos(lat) :",math.cos(math.radians(lat)))
if(event["origin"]["latitude"].get("uncertainty")):
laterr=float(event["origin"]["latitude"]["uncertainty"]) * 111.19 # from deg to km
print ("lat. uncert: ", laterr)
else:
print("NULL")
laterr=0.0
lon = float(event["origin"]["longitude"]["value"])
print(" long. ", lon)
if(event["origin"]["longitude"].get("uncertainty")):
lonerr=float(event["origin"]["longitude"]["uncertainty"]) \
* 111.19 * math.cos(math.radians(lat)) # from deg to km
print ( "long. uncert: ", lonerr)
else:
print("NULL")
lonerr=0.0
print("depth type: ",event["origin"]["depthType"])
if(event["origin"]["depth"].get("value")):
dep = float(event["origin"]["depth"]["value"]) /1000.0 # da metri a km
print("depth : ", dep)
else:
print("NULL")
dep = 10.0
if(event["origin"]["depth"].get("uncertainty")):
deperr = float(event["origin"]["depth"]["uncertainty"]) / 1000.0 # da metri a km
print ("depth uncert: ", deperr)
else:
print("NULL")
deperr = 5.0
if(event["origin"]["originUncertainty"].get("horizontalUncertainty")):
errh = float(event["origin"]["originUncertainty"]["horizontalUncertainty"]) / 1000.0 # da metri a km
print("horiz. uncert: ", errh)
else:
print("NULL")
errh = 5.0
if(event["origin"]["quality"].get("usedStationCount")):
usedsta = float(event["origin"]["quality"]["usedStationCount"])
print(" used stations: ", usedsta)
else:
print("NULL")
usedsta = 3.
if(event["origin"]["quality"].get("usedPhaseCount")):
usedpha = float(event["origin"]["quality"]["usedPhaseCount"])
print(" used phases: ", usedpha)
else:
print("NULL")
usedpha = 3.
if(event["origin"]["quality"].get("azimuthalGap")):
azimgap = float(event["origin"]["quality"]["azimuthalGap"])
print(" azim gap: ", azimgap )
else:
print("NULL")
azimgap = 350.
if(event["origin"]["quality"].get("minimumDistance")):
mindist = float(event["origin"]["quality"]["minimumDistance"]) * 111.19 # from deg to km
print(" min dist: ", mindist)
else:
print("NULL")
mindist = 50.
if(event["origin"]["quality"].get("maximumDistance")):
maxdist = float(event["origin"]["quality"]["maximumDistance"]) * 111.19 # from deg to km
print(" max dist: ", maxdist)
else:
print("NULL")
maxdist = 1000.
if(event["origin"]["quality"].get("standardError")):
stderr = float(event["origin"]["quality"]["standardError"])
print(" standard error: ", stderr)
else:
print("NULL")
stderr = 3.0
if(event["origin"]["creationInfo"].get("version")):
version = event["origin"]["creationInfo"]["version"]
print(" version: ", version)
else:
print("NULL")
version = "000"
print(event["description"]["text"])
if( event.get("magnitude") and event["magnitude"]["mag"].get("value")):
mag = float(event["magnitude"]["mag"]["value"])
magtype = event["magnitude"]["type"]
print("magnitude: ", magtype," ",mag)
else:
print("NO MAG !")
mag = -9.9
magtype = "NOMAG"
# applico i CRITERI ORIGINALI di Hypoinverse2000, con gli scores di Amato-Mele (2008)
#
# The first quality rating is based on errors and goodness-of-fit:
#
# 1.5 A. RMS < 0.15 sec and ERH <1.0 km and ERZ <2.0 km
# 0.5 B. RMS < 0.30 sec and ERH <2.5 km and ERZ <5.0 km
# -0.5 C. RMS < 0.50 sec and ERH <5.0 km
# -1.5 D. Worse than above
#
# The second quality rating is based on station geometry:
#
# 3.0 A. NWR > 6 and MAXGAP < 90 and either DMIN < DEPTH or DMIN <5.0
# 1.0 B. NWR > 6 and MAXGAP < 135 and either DMIN < 2*DEPTH or DMIN <10
# -1.0 C. NWR > 6 and MAXGAP < 180 and DMIN <100
# -3.0 D. Worse than above
# primo criterio:
if( stderr < 0.15 and errh < 1.0 and deperr < 2.0 ):
score1 = 1.5
elif( stderr < 0.30 and errh < 2.5 and deperr < 5.0 ):
score1 = 0.5
elif( stderr < 0.50 and errh < 5.0 and deperr < 5.0 ):
score1 = -0.5
else:
score1 = -1.5
# secondo criterio:
# Non ho a disposizione i pesi delle fasi, perciò adotto un criterio
# diverso: il criterio ricchiederebbe che la somma dei pesi delle fasi usate sia 6.
# Supponendo un valore medio dei pesi attribuito alle fasi di 0.75, richiedo che il numero
# totale di fasi sia almeno 8.
if( usedpha > 8 and azimgap < 90. and ( mindist < dep or mindist < 5.0)):
score2 = 3.0
elif( usedpha > 8 and azimgap < 135. and ( mindist < 2*dep or mindist < 10.0)):
score2 = 1.0
elif( usedpha > 8 and azimgap < 180. and mindist < 100.0):
score2 = -1.0
else:
score2 = -3.0
score = score1 + score2
print(" score1 score2 scoretot ", score1," ", score2," ", score)
# i codici GMT sono t, s, n, h, g, a (per i limiti di depth: 0 -> 5, 10, 30, 60, 100, > 100)
if ( dep <= 5):
dsig = "t"
elif( dep <= 10):
dsig = "s"
elif( dep <= 30):
dsig = "n"
elif( dep <= 60):
dsig = "h"
elif( dep <= 100):
dsig = "g"
else:
dsig = "a"
radius = mag*mag/40.0
print( "%.5f" % round(lon,5)," ",\
"%.5f" % round(lat,5)," ",\
"%.1f" % round(score,1)," ",\
"%.5f" % round(radius,5)," ",\
dsig, file=f)
# ....................=============================================
else:
print("version: ",version)
# ======================================================================
else:
print("version: NULL")
version = "000"
i=i+1
print("--------------------------------------",i)
f.close()
if __name__ == "__main__":
main()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment