Skip to content

Commit fb9939d

Browse files
strawpantsstrawpants
committed
move static gravity fields, degree 1 and 2 corrrections from geoslurp to this module
1 parent e540b20 commit fb9939d

File tree

5 files changed

+703
-0
lines changed

5 files changed

+703
-0
lines changed

‎pyproject.toml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,4 +39,8 @@ shlib = "shxarray.shlib:SHComputeBackend"
3939
"Homepage" = "https://github.com/ITC-Water-Resources/shxarray"
4040
"Bug Tracker" = "https://github.com/ITC-Water-Resources/shxarray/issues"
4141

42+
[project.entry-points."geoslurp.dsetfactories"]
43+
deg1n2corr = "shxarray.geoslurp.deg1n2:getDeg1n2corrDsets"
4244

45+
[project.entry-points."geoslurp.dsets"]
46+
icgemstatic = "shxarray.geoslurp.icgemdset:ICGEMstatic"

‎src/shxarray/geoslurp/deg1n2.py‎

Lines changed: 365 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,365 @@
1+
# This file is part of the shxarray software which is licensed
2+
# under the Apache License version 2.0 (see the LICENSE file in the main repository)
3+
# Copyright Roelof Rietbroek ([email protected]), 2024
4+
#
5+
6+
7+
from geoslurp.dataset import DataSet
8+
from geoslurp.datapull.ftp import Uri as ftp
9+
from geoslurp.datapull.http import Uri as http
10+
from datetime import datetime,timedelta
11+
from shxarray.core.time import decyear2dt
12+
import os
13+
import re
14+
from shxarray.geoslurp.gravity import GravitySHTBase,GravitySHinDBTBase,Trig,JSONSHArchive
15+
import xarray as xr
16+
17+
schema="shxarray"
18+
19+
# class geocenter_Rietbroeketal2016upd(DataSet):
20+
# fout="Geocenter_dec2017.tgz"
21+
# sqrt3timesRE=11047256.4063275
22+
# schema=schema
23+
# table=type("geocenter_Rietbroeketal2016updTable", (GravitySHTBase,), {})
24+
# def __init__(self,dbconn):
25+
# super().__init__(dbconn)
26+
# # super().__init__(direc=direc,uri='https://wobbly.earth/data/Geocenter_dec2017.tgz',order=['c10','c11','s11'],lastupdate=datetime(2018,10,16))
27+
28+
# def pull(self):
29+
# """Pulls the geocenter ascii files in the cache"""
30+
31+
# uri=http("https://wobbly.earth/data/"+self.fout,lastmod=datetime(2018,10,16)).download(self.cacheDir(),check=True)
32+
33+
34+
# def register(self):
35+
# self.truncateTable()
36+
# #set general settings
37+
# self._dbinvent.data={"citation":"Rietbroek, R., Brunnabend, S.-E., Kusche, J., Schröter, J., Dahle, C., 2016. " \
38+
# "Revisiting the Contemporary Sea Level Budget on Global and Regional Scales. " \
39+
# "Proceedings of the National Academy of Sciences 201519132. " \
40+
# "https://doi.org/10.1073/pnas.1519132113"}
41+
42+
# with tarfile.open(os.path.join(self.cacheDir(),self.fout),'r:gz') as tf:
43+
44+
# metacomb=[]
45+
46+
# files=['Geocenter/GeocentCM-CF_Green.txt',
47+
# 'Geocenter/GeocentCM-CF_Antarctica.txt',
48+
# 'Geocenter/GeocentCM-CF_Hydrology.txt',
49+
# 'Geocenter/GeocentCM-CF_LandGlaciers.txt',
50+
# 'Geocenter/GeocentCM-CF_GIA.txt',
51+
# 'Geocenter/GeocentCM-CF_TotSurfload.txt']
52+
53+
# order=[(1,1,Trig.c),(1,1,Trig.s),(1,0,Trig.c)]
54+
# lastupdate=datetime.now()
55+
# for file in files:
56+
# #get files
57+
# with tf.extractfile(file) as fid:
58+
# for ln in fid:
59+
# shar=JSONSHArchive(1)
60+
# lnspl=ln.decode('utf-8').split()
61+
# tcent=decyear2dt(float(lnspl[0]))
62+
# tstart,tend=dt2monthlyinterval(tcent)
63+
64+
# meta={"type":file.split('_')[-1][:-4],"time":tcent,"tstart":tstart,"tend":tend,"lastupdate":lastupdate,"nmax":1,"omax":1,"origin":"CF","format":"JSONB","uri":"self:data","gm":0.3986004415e+15,"re":0.6378136460e+07
65+
# }
66+
67+
# for el,val in zip(order,lnspl[1:4]):
68+
# # import pdb;pdb.set_trace()
69+
# shar["cnm"][shar.idx(el)]=float(val)/self.sqrt3timesRE
70+
71+
# #also add sigmas
72+
# for el,val in zip(order,lnspl[4:7]):
73+
# shar["sigcnm"][shar.idx(el)]=float(val)/self.sqrt3timesRE
74+
# meta["data"]=shar.dict
75+
# self.addEntry(meta)
76+
# self.updateInvent()
77+
78+
79+
def parseGSMDate(dtstr):
80+
"""Parse datestr as found in GSM files (yyyymmdd.00000)"""
81+
return datetime(int(dtstr[0:4]),int(dtstr[4:6]),int(dtstr[6:8]))
82+
83+
class geocenter_GRCRL06_TN13(DataSet):
84+
schema=schema
85+
rooturl="http://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-docs/grace/open/docs/"
86+
def __init__(self,dbconn):
87+
self.table=type(self.__class__.__name__.lower().replace('-',"_")+"Table", (GravitySHinDBTBase,), {})
88+
super().__init__(dbconn)
89+
90+
def pull(self):
91+
"""Pulls the geocenter ascii files in the cache"""
92+
uri=http(self.rooturl+self.fout,lastmod=datetime(2019,12,1)).download(self.cacheDir(),check=True)
93+
94+
95+
def register(self):
96+
self.truncateTable()
97+
#set general settings
98+
self._dbinvent.data={"citation":"GRACE technical note 13"}
99+
lastupdate=datetime.now()
100+
with open(os.path.join(self.cacheDir(),self.fout),'r') as fid:
101+
#skip header
102+
for ln in fid:
103+
if re.search("end of header",ln):
104+
break
105+
106+
nv=[]
107+
mv=[]
108+
tv=[]
109+
cnmv=[]
110+
sigcnmv=[]
111+
#loop over entry lines
112+
for ln in fid:
113+
114+
tp,n,m,cnm,snm,sigcnm,sigsnm,ts,te=ln.split()
115+
116+
#Append cosine coefficients
117+
n=int(n)
118+
m=int(m)
119+
nv.append(n)
120+
mv.append(m)
121+
tv.append(0)
122+
cnmv.append(float(cnm))
123+
sigcnmv.append(float(sigcnm))
124+
125+
#append sine coefficients
126+
if m > 0:
127+
nv.append(n)
128+
mv.append(m)
129+
tv.append(1)
130+
cnmv.append(float(snm))
131+
sigcnmv.append(float(sigsnm))
132+
133+
if m == 1:
134+
#register the accumulated entry
135+
tstart=parseGSMDate(ts)
136+
tend=parseGSMDate(te)
137+
#get the central time
138+
tcent=tstart+(tend-tstart)/2
139+
#snap the central epoch to the 15th of the month of the central time
140+
# tcent=datetime(tstart.year,tstart.month,15)
141+
142+
meta={"type":"GSM","time":tcent,"tstart":tstart,"tend":tend,"lastupdate":lastupdate,"nmax":1,"omax":1,"origin":"CF","format":"JSONB","gm":0.3986004415e+15,"re":0.6378136460e+07}
143+
meta["data"]=xr.Dataset(data_vars=dict(cnm=(["shg"],cnmv),sigcnm=(["shg"],sigcnmv)),coords=dict(n=(["shg"],nv),m=(["shg"],mv),t=(["shg"],tv)))
144+
145+
self.addEntry(meta)
146+
nv=[]
147+
mv=[]
148+
tv=[]
149+
cnmv=[]
150+
sigcnmv=[]
151+
152+
self.updateInvent()
153+
154+
155+
156+
157+
158+
159+
class GeocenterRIESCFCM(DataSet):
160+
fout30="GCN_L1_L2_30d_CF-CM.txt"
161+
fout60="GCN_L1_L2_60d_CF-CM.txt"
162+
#note also embed mm to meter conversion in here (e3)
163+
sqrt3timesRE=11047256.23312e3
164+
schema=schema
165+
table=type("geocenter_ries_cfcmTable", (GravitySHTBase,), {})
166+
def __init__(self,dbconn):
167+
super().__init__(dbconn)
168+
# super().__init__(direc=direc,uri='https://wobbly.earth/data/Geocenter_dec2017.tgz',order=['c10','c11','s11'],lastupdate=datetime(2018,10,16))
169+
170+
def pull(self):
171+
"""Pulls the geocenter ascii files in the cache"""
172+
173+
uri=http("http://download.csr.utexas.edu/pub/slr/geocenter/"+self.fout30).download(self.cacheDir())
174+
uri=http("http://download.csr.utexas.edu/pub/slr/geocenter/"+self.fout60).download(self.cacheDir())
175+
176+
177+
def register(self):
178+
self.truncateTable()
179+
#set general settings
180+
self._dbinvent.data={"citation":"Ries, J.C., 2016. Reconciling estimates of annual geocenter motion from space geodesy, in: Proceedings of the 20th International Workshop on Laser Ranging, Potsdam, Germany. pp. 10–14."}
181+
182+
183+
self.extractSLR(os.path.join(self.cacheDir(),self.fout30))
184+
self.extractSLR(os.path.join(self.cacheDir(),self.fout60))
185+
self.updateInvent()
186+
187+
def extractSLR(self,filen):
188+
lastupdate=datetime.now()
189+
order=[(1,1,Trig.c),(1,1,Trig.s),(1,0,Trig.c)]
190+
if re.search('30d',filen):
191+
dt=timedelta(days=15)
192+
else:
193+
dt=timedelta(days=30)
194+
195+
with open(filen,'r') as fid:
196+
#skip header
197+
fid.readline()
198+
for ln in fid:
199+
shar=JSONSHArchive(1)
200+
lnspl=ln.split()
201+
#note little hack as the 60d file version has an empty line at the back
202+
if len(lnspl) == 0:
203+
break
204+
tcent=decyear2dt(float(lnspl[0]))
205+
tstart=tcent-dt
206+
tend=tcent+dt
207+
208+
meta={"type":"GSM"+os.path.basename(filen)[10:13],"time":tcent,"tstart":tstart,"tend":tend,"lastupdate":lastupdate,"nmax":1,"omax":1,"origin":"CF","format":"JSONB","uri":"self:data","gm":0.3986004415e+15,"re":0.6378136460e+07
209+
}
210+
211+
for el,val in zip(order,lnspl[1:4]):
212+
shar["cnm"][shar.idx(el)]=float(val)/self.sqrt3timesRE
213+
214+
#also add sigmas
215+
for el,val in zip(order,lnspl[4:7]):
216+
shar["sigcnm"][shar.idx(el)]=float(val)/self.sqrt3timesRE
217+
meta["data"]=shar.dict
218+
self.addEntry(meta)
219+
220+
221+
class TN14SLRGSFC(DataSet):
222+
schema=schema
223+
rooturl="ftp://isdcftp.gfz-potsdam.de/grace-fo/DOCUMENTS/TECHNICAL_NOTES/"
224+
fout="TN-14_C30_C20_SLR_GSFC.txt"
225+
def __init__(self,dbconn):
226+
self.table=type(self.__class__.__name__.lower().replace('-',"_")+"Table", (GravitySHinDBTBase,), {})
227+
super().__init__(dbconn)
228+
229+
def pull(self):
230+
"""Pulls the C20 Technical note 14"""
231+
uri=ftp(self.rooturl+self.fout,lastmod=datetime(2019,12,1)).download(self.cacheDir(),check=True)
232+
233+
234+
def register(self):
235+
self.truncateTable()
236+
#set general settings
237+
self._dbinvent.data={"citation":"GRACE technical note 14"}
238+
lastupdate=datetime.now()
239+
mjd00=datetime(1858,11,17)
240+
nmax=2
241+
omax=0
242+
with open(os.path.join(self.cacheDir(),self.fout),'r') as fid:
243+
#skip header
244+
for ln in fid:
245+
if re.search("^Product:",ln):
246+
break
247+
248+
#loop over entry lines
249+
for ln in fid:
250+
251+
mjd0,decy0,c20,dc20,sigc20,c30,dc30,sigc30,mjd1,decyr1=ln.split()
252+
253+
nv=[]
254+
mv=[]
255+
tv=[]
256+
cnmv=[]
257+
dcnmv=[]
258+
sigcnmv=[]
259+
260+
#Append c20 coefficients
261+
nv.append(2)
262+
mv.append(0)
263+
tv.append(0)
264+
cnmv.append(float(c20))
265+
dcnmv.append(float(dc20)*1e-10)
266+
sigcnmv.append(float(sigc20))
267+
if c30 != "NaN":
268+
nmax=3
269+
nv.append(3)
270+
mv.append(0)
271+
tv.append(0)
272+
cnmv.append(float(c30))
273+
dcnmv.append(float(dc30)*1e-10)
274+
sigcnmv.append(float(sigc30))
275+
276+
#register the accumulated entry
277+
tstart=mjd00+timedelta(days=float(mjd0))
278+
tend=mjd00+timedelta(days=float(mjd1))
279+
tcent=tstart+(tend-tstart)/2
280+
281+
meta={"type":"GSM","time":tcent,"tstart":tstart,"tend":tend,"lastupdate":lastupdate,"nmax":nmax,"omax":omax,"format":"JSONB","gm":0.3986004415e+15,"re":0.6378136460e+07}
282+
meta["data"]=xr.Dataset(data_vars=dict(cnm=(["shg"],cnmv),dcnm=(["shg"],dcnmv),sigcnm=(["shg"],sigcnmv)),coords=dict(n=(["shg"],nv),m=(["shg"],mv),t=(["shg"],tv)))
283+
284+
self.addEntry(meta)
285+
self.updateInvent()
286+
287+
def getDeg1n2corrDsets(conf):
288+
out=[]
289+
for center in ["CSR", "GFZ", "JPL"]:
290+
out.append(type("geocenter_"+center+"RL06_TN13",(geocenter_GRCRL06_TN13,),{"fout":"TN-13_GEOC_"+center+"_RL06.txt"}))
291+
out.append(GeocenterRIESCFCM)
292+
out.append(TN14SLRGSFC)
293+
return out
294+
295+
# class Sun2017Comb(LowdegreeSource):
296+
# def __init__(self,direc,uri=None):
297+
# if not uri:
298+
# uri='https://d1rkab7tlqy5f1.cloudfront.net/CiTG/Over%20faculteit/Afdelingen/' \
299+
# 'Geoscience%20%26%20Remote%20sensing/Research/Gravity/models%20and%20data/3_Deg1_C20_CMB.txt'
300+
# lastupdate=datetime(2018,1,1)
301+
# super().__init__(direc=direc,uri=uri,lastupdate=lastupdate)
302+
# def extract(self):
303+
# singlemeta=self.meta
304+
305+
# citation="Reference: Yu Sun, Pavel Ditmar, Riccardo Riva (2017), Statistically optimal" \
306+
# " estimation of degree-1 and C20 coefficients based on GRACE data and an " \
307+
# "ocean bottom pressure model Geophysical Journal International, 210(3), " \
308+
# "1305-1322, 2017. doi:10.1093/gji/ggx241."
309+
310+
311+
# valorder=['c10','c11','s11']
312+
# covorder=[('c10','c10'),('c10','c11'),('c10','s11'),('c10','c20'),('c11','c11'),('c11','s11'),('c11','c20'),('s11','s11'),('s11','c20'),('c20','c20')]
313+
314+
# singlemeta['data']["citation"]=citation
315+
# with open(os.path.join(self.direc,self.fout),'rt') as fid:
316+
# dataregex=re.compile('^ +[0-9]')
317+
# time=[]
318+
# for ln in fid:
319+
# if dataregex.match(ln):
320+
# lnspl=ln.split()
321+
# time.append(decyear2dt(float(lnspl[0])))
322+
# d12=[0]*4
323+
# for el,val in zip(valorder,lnspl[1:4]):
324+
# d12[self.dindex(el)]=val
325+
# singlemeta['data']['d12'].append(d12)
326+
327+
328+
# covUpper=[0]*10
329+
# for (el1,el2),val in zip(covorder,lnspl[5:]):
330+
# covUpper[self.covindex(el1,el2)]=val
331+
# singlemeta['data']['covUpper'].append(covUpper)
332+
333+
# singlemeta["data"]["time"]=[dt.isoformat() for dt in time]
334+
# singlemeta['tstart']=min(time)
335+
# singlemeta['tend']=max(time)
336+
# singlemeta['lastupdate']=datetime.now()
337+
338+
# return [singlemeta]
339+
340+
# class Sun2017Comb_GIArestored(Sun2017Comb):
341+
# def __init__(self,direc):
342+
# uri="https://d1rkab7tlqy5f1.cloudfront.net/CiTG/Over%20faculteit/Afdelingen/" \
343+
# "Geoscience%20%26%20Remote%20sensing/Research/Gravity/models%20and%20data/" \
344+
# "4_Deg1_C20_CMB_GIA_restored.txt"
345+
# super().__init__(direc,uri)
346+
347+
# class GeocTable(GeocTBase):
348+
# """Defines the Geocenter motion table"""
349+
# __tablename__='deg1n2'
350+
# id=Column(Integer,primary_key=True)
351+
# name=Column(String,unique=True)
352+
# lastupdate=Column(TIMESTAMP)
353+
# tstart=Column(TIMESTAMP)
354+
# tend=Column(TIMESTAMP)
355+
# origin=Column(String)
356+
# data=Column(JSONB)
357+
358+
# "ftp://ftp.csr.utexas.edu/pub/slr/geocenter/"
359+
# # {'name':'Rietbroeketal2016updated','uri':'https://wobbly.earth/data/Geocenter_dec2017.tgz','lastupdate':datetime(2018,10,16)},
360+
# # {'name':'SwensonWahr2008','uri':'ftp://podaac.jpl.nasa.gov/allData/tellus/L2/degree_1/deg1_coef.txt','lastupdate':datetime(2018,10,16)},
361+
# # {'name':'Sun2017Comb','uri':'https://d1rkab7tlqy5f1.cloudfront.net/CiTG/Over%20faculteit/Afdelingen/Geoscience%20%26%20Remote%20sensing/Research/Gravity/models%20and%20data/3_Deg1_C20_CMB.txt'
362+
# # ,'lastupdate':datetime(2018,10,16)},
363+
# ]
364+
365+

0 commit comments

Comments
 (0)