From ebe557cc3444c3df86d839abc3016c5adeac56c1 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 14 Apr 2023 11:09:39 -0400 Subject: [PATCH 1/7] convert solid_point to return an array instead of writing to text file --- src/pysolid/point.py | 20 ++------------- src/pysolid/solid.for | 57 +++++++++++++++++-------------------------- 2 files changed, 24 insertions(+), 53 deletions(-) diff --git a/src/pysolid/point.py b/src/pysolid/point.py index 07d109a..37c1924 100644 --- a/src/pysolid/point.py +++ b/src/pysolid/point.py @@ -169,23 +169,10 @@ def calc_solid_earth_tides_point_per_day(lat, lon, date_str, step_sec=60): msg += '\n Check instruction at: https://github.com/insarlab/PySolid.' raise ImportError(msg) - ## calc solid Earth tides and write to text file - txt_file = os.path.abspath('solid.txt') - if os.path.isfile(txt_file): - os.remove(txt_file) - + # calc solid Earth tides # Run twice to circumvent fortran bug which cuts off last file in loop - Simran, Jun 2020 t = dt.datetime.strptime(date_str, '%Y%m%d') - for _ in range(2): - solid_point(lat, lon, t.year, t.month, t.day, step_sec) - - ## read data from text file - num_row = int(24 * 60 * 60 / step_sec) - fc = np.loadtxt(txt_file, - dtype=float, - delimiter=',', - skiprows=0, - max_rows=num_row) + fc = solid_point(lat, lon, t.year, t.month, t.day, step_sec) tide_e = fc[:, 1].flatten() tide_n = fc[:, 2].flatten() @@ -195,9 +182,6 @@ def calc_solid_earth_tides_point_per_day(lat, lon, date_str, step_sec=60): dt_out = [t + dt.timedelta(seconds=sec) for sec in secs] dt_out = np.array(dt_out) - # remove the temporary text file - os.remove(txt_file) - return dt_out, tide_e, tide_n, tide_u diff --git a/src/pysolid/solid.for b/src/pysolid/solid.for index 4cf31bd..6c2dd6b 100644 --- a/src/pysolid/solid.for +++ b/src/pysolid/solid.for @@ -8,7 +8,6 @@ *** dependent subroutines can be download from IERS conventions website as: *** wget -r -l1 --no-parent -R "index.html*" -nH --cut-dirs=3 https://iers-conventions.obspm.fr/content/chapter7/software/dehanttideinel *** Z. Yunjun and S. Sangha, Sep 2020: modify solid() to solid_point/grid() as subroutines. - subroutine solid_grid(iyr,imo,idy,ihh,imm,iss, * glad0,steplat,nlat, * glod0,steplon,nlon) @@ -174,13 +173,14 @@ 99 end *----------------------------------------------------------------------- - subroutine solid_point(glad,glod,iyr,imo,idy,step_sec) + subroutine solid_point(glad,glod,iyr,imo,idy,step_sec,output) *** calculate SET at given location for one day with step_sec seconds resolution *** Arguments: glad/glod - float, latitude/longitude in deg *** iyr/imo/idy - int, start date/time in UTC *** step_sec - int, time step in seconds -*** Returns: seconds, SET_east, SET_north, SET_up +*** Returns: output +*** each row: seconds, SET_east, SET_north, SET_up implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -188,16 +188,13 @@ integer iyr,imo,idy integer nloop, step_sec double precision tdel2 + real(8), intent(out), dimension(60*60*24/step_sec,4) :: output !*** leap second table limit flag logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 - -*** open output file - - lout=1 - open(lout,file='solid.txt',form='formatted',status='unknown') - write(lout,'(a)') '# program solid -- UTC version -- 2018jun01' +! f2py intent(in) glad,glod,iyr,imo,idy,nloop +! f2py intent(out) output *** constants @@ -213,34 +210,30 @@ *** check inputs section if(glad.lt.-90.d0.or.glad.gt.90.d0) then - write(lout,'(a,G0.9)') 'ERROR: lat NOT in [-90,+90]:',glad - go to 98 + print *, 'ERROR: lat NOT in [-90,+90]:',glad + return endif if(glod.lt.-360.d0.or.glod.gt.360.d0) then - write(lout,'(a,G0.9)') 'ERROR: lon NOT in [-360,+360]:',glod - go to 98 + print *, 'ERROR: lon NOT in [-360,+360]:',glod + return endif if(iyr.lt.1901.or.iyr.gt.2099) then - write(lout,'(a,i5)') 'ERROR: year NOT in [1901-2099]:',iyr - go to 98 + print *, 'ERROR: year NOT in [1901-2099]:',iyr + return endif if(imo.lt.1.or.imo.gt.12) then - write(lout,'(a,i3)') 'ERROR: month NOT in [1-12]:',imo - go to 98 + print *, 'ERROR: month NOT in [1-12]:',imo + return endif if(idy.lt.1.or.idy.gt.31) then - write(lout,'(a,i3)') 'ERROR: day NOT in [1-31]:',idy - go to 98 + print *, 'ERROR: day NOT in [1-31]:',idy + return endif -*** output header - - write(lout,'(a,i5,2i3)') '# year, month, day =',iyr,imo,idy - write(lout,'(a,2f15.9)') '# lat, lon =',glad,glod *** position of observing point (positive East) @@ -268,11 +261,9 @@ *** loop over time - !*** tdel2=1.d0/60.d0/24.d0 - !*** do iloop=0,60*24 nloop=60*60*24/step_sec tdel2=1.d0/DFLOAT(nloop) - do iloop=0,nloop + do iloop=1,nloop !*** false means flag not raised !*** mjd/fmjd in UTC !*** mjd/fmjd in UTC @@ -292,9 +283,9 @@ call mjdciv(mjd,fmjd,iyr,imo,idy,ihr,imn,sec) - !*** write output to file tsec=ihr*3600.d0+imn*60.d0+sec - write(lout,'((f8.1,:,", "),*(f10.6,:,", "))') tsec,vt,ut,wt + + output(iloop, :) = [tsec,vt,ut,wt] !*** update fmjd for the next round fmjd=fmjd+tdel2 @@ -305,16 +296,12 @@ *** end of processing and flag of leap second if(lflag) then - write(*,'(a)') 'Mild Warning -- time crossed leap second table' - write(*,'(a)') ' boundaries. Boundary edge value used instead' + print *, 'Mild Warning -- time crossed leap second table' + print *, ' boundaries. Boundary edge value used instead' endif - close(lout) - - go to 99 - 98 write(*,'(a)') 'Check arguments.' return - 99 end + end *@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ subroutine detide(xsta,mjd,fmjd,xsun,xmon,dxtide,lflag) From 29c968c97f349f414b6bfa6af9c51254c8ed9b08 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 14 Apr 2023 13:56:36 -0400 Subject: [PATCH 2/7] change solid grid fortran to return arrays. fixes #2, fixes #50 try extra line length fix space/multiline issue of long subroutine --- setup.py | 4 ++- src/pysolid/grid.py | 37 ++++++-------------- src/pysolid/solid.for | 81 +++++++++++++++++++------------------------ 3 files changed, 48 insertions(+), 74 deletions(-) diff --git a/setup.py b/setup.py index 8b49253..c56034a 100644 --- a/setup.py +++ b/setup.py @@ -58,6 +58,8 @@ def readme(): ## fortran extensions to build with numpy.f2py ext_modules=[ - Extension(name="pysolid.solid", sources=["src/pysolid/solid.for"]), + Extension( + name="pysolid.solid", sources=["src/pysolid/solid.for"], + ) ], ) diff --git a/src/pysolid/grid.py b/src/pysolid/grid.py index 3920943..f2eba23 100644 --- a/src/pysolid/grid.py +++ b/src/pysolid/grid.py @@ -73,35 +73,18 @@ def calc_solid_earth_tides_grid(dt_obj, atr, step_size=1e3, display=False, verbo vprint('SOLID : shape: {s}, step size: {la:.4f} by {lo:.4f} deg'.format( s=(length, width), la=lat_step, lo=lon_step)) - ## calc solid Earth tides and write to text file - txt_file = os.path.abspath('solid.txt') - if os.path.isfile(txt_file): - os.remove(txt_file) + ## calc solid Earth tides + # Run twice to circumvent fortran bug which cuts off last file in loop - Simran, Jun 2020 + fc = solid_grid(dt_obj.year, dt_obj.month, dt_obj.day, + dt_obj.hour, dt_obj.minute, dt_obj.second, + lat0, lat_step, length, + lon0, lon_step, width) - vprint('SOLID : calculating / writing data to txt file: {}'.format(txt_file)) + tide_e = fc[:, :, 0].reshape(length, width) + tide_n = fc[:, :, 1].reshape(length, width) + tide_u = fc[:, :, 2].reshape(length, width) - # Run twice to circumvent fortran bug which cuts off last file in loop - Simran, Jun 2020 - for _ in range(2): - solid_grid(dt_obj.year, dt_obj.month, dt_obj.day, - dt_obj.hour, dt_obj.minute, dt_obj.second, - lat0, lat_step, length-1, - lon0, lon_step, width-1) - - ## read data from text file - vprint('PYSOLID: read data from text file: {}'.format(txt_file)) - grid_size = int(length * width) - fc = np.loadtxt(txt_file, - dtype=float, - usecols=(2,3,4), - delimiter=',', - skiprows=0, - max_rows=grid_size+100)[:grid_size] - tide_e = fc[:, 0].reshape(length, width) - tide_n = fc[:, 1].reshape(length, width) - tide_u = fc[:, 2].reshape(length, width) - - # remove the temporary text file - os.remove(txt_file) + np.save('fc_grid_arr.npy', fc) # resample to the input size if num_step > 1: diff --git a/src/pysolid/solid.for b/src/pysolid/solid.for index 6c2dd6b..8c0bd3f 100644 --- a/src/pysolid/solid.for +++ b/src/pysolid/solid.for @@ -9,14 +9,14 @@ *** wget -r -l1 --no-parent -R "index.html*" -nH --cut-dirs=3 https://iers-conventions.obspm.fr/content/chapter7/software/dehanttideinel *** Z. Yunjun and S. Sangha, Sep 2020: modify solid() to solid_point/grid() as subroutines. subroutine solid_grid(iyr,imo,idy,ihh,imm,iss, - * glad0,steplat,nlat, - * glod0,steplon,nlon) - + * glad0,steplat,nlat,glod0,steplon,nlon,output) + *** calculate solid earth tides (SET) for one spatial grid given the date/time *** Arguments: iyr/imo/idy/ihh/imm/iss - int, date/time for YYYY/MM/DD/HH/MM/SS *** glad0/glad1/steplat - float, north(Y_FIRST)/south/step(negative) in deg *** glod0/glod1/steplon - float, west(X_FIRST) /east /step(positive) in deg -*** Returns: latitude, longitude, SET_east, SET_north, SET_up +*** Returns: output: 3D array with shape (nlat, nlon, 5) +*** each depth slice: [SET_east, SET_north, SET_up] implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -24,16 +24,14 @@ integer nlat,nlon double precision glad0,steplat double precision glod0,steplon + real(8), intent(out), dimension(nlat,nlon,3) :: output !***^ leap second table limit flag logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 - -*** open output file - - lout=1 - open(lout,file='solid.txt',form='formatted',status='unknown') - write(lout,'(a)') '# program solid -- UTC version -- 2018jun01' +! f2py intent(in) iyr,imo,idy,ihh,imm,iss,glad0,steplat,nlat,glod0,steplon,nlon +! f2py dimension(nlat,nlon,3) output +! f2py intent(out) output *** constants @@ -49,43 +47,43 @@ *** input section if(iyr.lt.1901.or.iyr.gt.2099) then - write(lout,'(a,i5)') 'ERROR: year NOT in [1901-2099]:',iyr - go to 98 + print *, 'ERROR: year NOT in [1901-2099]:',iyr + return endif if(imo.lt.1.or.imo.gt.12) then - write(lout,'(a,i3)') 'ERROR: month NOT in [1-12]:',imo - go to 98 + print *, 'ERROR: month NOT in [1-12]:',imo + return endif if(idy.lt.1.or.idy.gt.31) then - write(lout,'(a,i3)') 'ERROR: day NOT in [1-31]:',idy - go to 98 + print *, 'ERROR: day NOT in [1-31]:',idy + return endif if(ihh.lt.0.or.ihh.gt.23) then - write(lout,'(a,i3)') 'ERROR: hour NOT in [0-23]:',ihh - go to 98 + print *, 'ERROR: hour NOT in [0-23]:',ihh + return endif if(imm.lt.0.or.imm.gt.59) then - write(lout,'(a,i3)') 'ERROR: minute NOT in [0-59]:',imm - go to 98 + print *, 'ERROR: minute NOT in [0-59]:',imm + return endif if(iss.lt.0.or.iss.gt.59) then - write(lout,'(a,i3)') 'ERROR: second NOT in [0-59]:',iss - go to 98 + print *, 'ERROR: second NOT in [0-59]:',iss + return endif if(glad0.lt.-90.d0.or.glad0.gt.90.d0) then - write(lout,'(a,G0.9)') 'ERROR: lat0 NOT in [-90,+90]:',glad0 - go to 98 + print *, 'ERROR: lat0 NOT in [-90,+90]:',glad0 + return endif if(glod0.lt.-360.d0.or.glod0.gt.360.d0) then - write(lout,'(a,G0.9)') 'ERROR: lon0 NOT in [-360,+360]',glod0 - go to 98 + print *, 'ERROR: lon0 NOT in [-360,+360]',glod0 + return endif *** output header @@ -93,19 +91,13 @@ glad1=glad0+nlat*steplat glod1=glod0+nlon*steplon - write(lout,'(a,i5,2i3)') '# year, month, day =',iyr,imo,idy - write(lout,'(a,3i3)') '# hour, minute, second =',ihh,imm,iss - write(lout,'(a,4f15.9)') '# S, N, W, E =',glad1,glad0,glod0,glod1 - write(lout,'(a,f15.9,i6)') '# step_lat, num_lat =',steplat,nlat - write(lout,'(a,f15.9,i6)') '# step_lon, num_lon =',steplon,nlon - *** loop over the grid - do ilat=0,nlat - do ilon=0,nlon + do ilat=1,nlat + do ilon=1,nlon - glad=glad0+ilat*steplat - glod=glod0+ilon*steplon + glad = glad0 + (ilat-1)*steplat + glod = glod0 + (ilon-1)*steplon *** position of observing point (positive East) @@ -153,7 +145,7 @@ * iyr,imo,idy,ihr,imn,sec-0.001d0) !*** write output to file - write(lout,'(*(G0.9,:,", "))') glad,glod,vt,ut,wt + output(ilat, ilon, :) = [vt, ut, wt] enddo enddo @@ -161,16 +153,13 @@ *** end of processing and flag for leap second if(lflag) then - write(*,'(a)') 'Mild Warning -- time crossed leap second table' - write(*,'(a)') ' boundaries. Boundary edge value used instead' + print *, 'Mild Warning -- time crossed leap second table' + print *, ' boundaries. Boundary edge value used instead' endif - close(lout) - go to 99 - 98 write(*,'(a)') 'Check arguments.' return - 99 end + end *----------------------------------------------------------------------- subroutine solid_point(glad,glod,iyr,imo,idy,step_sec,output) @@ -180,7 +169,7 @@ *** iyr/imo/idy - int, start date/time in UTC *** step_sec - int, time step in seconds *** Returns: output -*** each row: seconds, SET_east, SET_north, SET_up +*** each row: [seconds, SET_east, SET_north, SET_up] implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -193,7 +182,7 @@ logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 -! f2py intent(in) glad,glod,iyr,imo,idy,nloop +! f2py intent(in) glad,glod,iyr,imo,idy,step_sec ! f2py intent(out) output *** constants @@ -285,7 +274,7 @@ tsec=ihr*3600.d0+imn*60.d0+sec - output(iloop, :) = [tsec,vt,ut,wt] + output(iloop, :) = [tsec, vt, ut, wt] !*** update fmjd for the next round fmjd=fmjd+tdel2 From ce72f430cd96cd6723fff920e097d64a05682f12 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 14 Apr 2023 15:20:49 -0400 Subject: [PATCH 3/7] remove "run twice" comment now that we run once --- setup.py | 4 +--- src/pysolid/grid.py | 3 --- src/pysolid/point.py | 1 - 3 files changed, 1 insertion(+), 7 deletions(-) diff --git a/setup.py b/setup.py index c56034a..9d5ef78 100644 --- a/setup.py +++ b/setup.py @@ -58,8 +58,6 @@ def readme(): ## fortran extensions to build with numpy.f2py ext_modules=[ - Extension( - name="pysolid.solid", sources=["src/pysolid/solid.for"], - ) + Extension(name="pysolid.solid", sources=["src/pysolid/solid.for"]) ], ) diff --git a/src/pysolid/grid.py b/src/pysolid/grid.py index f2eba23..89ceffe 100644 --- a/src/pysolid/grid.py +++ b/src/pysolid/grid.py @@ -74,7 +74,6 @@ def calc_solid_earth_tides_grid(dt_obj, atr, step_size=1e3, display=False, verbo s=(length, width), la=lat_step, lo=lon_step)) ## calc solid Earth tides - # Run twice to circumvent fortran bug which cuts off last file in loop - Simran, Jun 2020 fc = solid_grid(dt_obj.year, dt_obj.month, dt_obj.day, dt_obj.hour, dt_obj.minute, dt_obj.second, lat0, lat_step, length, @@ -84,8 +83,6 @@ def calc_solid_earth_tides_grid(dt_obj, atr, step_size=1e3, display=False, verbo tide_n = fc[:, :, 1].reshape(length, width) tide_u = fc[:, :, 2].reshape(length, width) - np.save('fc_grid_arr.npy', fc) - # resample to the input size if num_step > 1: out_shape = (int(atr['LENGTH']), int(atr['WIDTH'])) diff --git a/src/pysolid/point.py b/src/pysolid/point.py index 37c1924..af87fe6 100644 --- a/src/pysolid/point.py +++ b/src/pysolid/point.py @@ -170,7 +170,6 @@ def calc_solid_earth_tides_point_per_day(lat, lon, date_str, step_sec=60): raise ImportError(msg) # calc solid Earth tides - # Run twice to circumvent fortran bug which cuts off last file in loop - Simran, Jun 2020 t = dt.datetime.strptime(date_str, '%Y%m%d') fc = solid_point(lat, lon, t.year, t.month, t.day, step_sec) From f6d8de2ab0bec721ed0476169e1d23a290c7f911 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 16 Apr 2023 09:00:27 -0400 Subject: [PATCH 4/7] fix bad merging of main branch --- src/pysolid/grid.py | 2 +- src/pysolid/point.py | 1 - src/pysolid/solid.for | 10 +++------- 3 files changed, 4 insertions(+), 9 deletions(-) diff --git a/src/pysolid/grid.py b/src/pysolid/grid.py index 2aa7941..89ceffe 100644 --- a/src/pysolid/grid.py +++ b/src/pysolid/grid.py @@ -11,8 +11,8 @@ # pysolid.calc_solid_earth_tides_grid() +import datetime as dt import os -import tempfile import numpy as np from skimage.transform import resize diff --git a/src/pysolid/point.py b/src/pysolid/point.py index bd2e457..af87fe6 100644 --- a/src/pysolid/point.py +++ b/src/pysolid/point.py @@ -14,7 +14,6 @@ import collections import datetime as dt import os -import tempfile import numpy as np from matplotlib import pyplot as plt, ticker, dates as mdates diff --git a/src/pysolid/solid.for b/src/pysolid/solid.for index e17cc69..7188ccc 100644 --- a/src/pysolid/solid.for +++ b/src/pysolid/solid.for @@ -10,17 +10,15 @@ *** Z. Yunjun and S. Sangha, Sep 2020: modify solid() to solid_point/grid() as subroutines. subroutine solid_grid(iyr,imo,idy,ihh,imm,iss, * glad0,steplat,nlat,glod0,steplon,nlon,output) - + *** calculate solid earth tides (SET) for one spatial grid given the date/time -*** Arguments: txt_file - string, output file name -*** iyr/imo/idy/ihh/imm/iss - int, date/time for YYYY/MM/DD/HH/MM/SS +*** Arguments: iyr/imo/idy/ihh/imm/iss - int, date/time for YYYY/MM/DD/HH/MM/SS *** glad0/glad1/steplat - float, north(Y_FIRST)/south/step(negative) in deg *** glod0/glod1/steplon - float, west(X_FIRST) /east /step(positive) in deg *** Returns: output: 3D array with shape (nlat, nlon, 5) *** each depth slice: [SET_east, SET_north, SET_up] implicit double precision(a-h,o-z) - character(len=*), intent(in) :: txt_file dimension rsun(3),rmoon(3),etide(3),xsta(3) integer iyr,imo,idy,ihh,imm,iss integer nlat,nlon @@ -167,15 +165,13 @@ subroutine solid_point(glad,glod,iyr,imo,idy,step_sec,output) *** calculate SET at given location for one day with step_sec seconds resolution -*** Arguments: txt_file - string, output file name -*** glad/glod - float, latitude/longitude in deg +*** Arguments: glad/glod - float, latitude/longitude in deg *** iyr/imo/idy - int, start date/time in UTC *** step_sec - int, time step in seconds *** Returns: output *** each row: [seconds, SET_east, SET_north, SET_up] implicit double precision(a-h,o-z) - character(len=*), intent(in) :: txt_file dimension rsun(3),rmoon(3),etide(3),xsta(3) double precision glad,glod integer iyr,imo,idy From 9548b9099931f8a5b309e107ff6a715e2a2fbe9b Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 16 Apr 2023 16:09:24 -0400 Subject: [PATCH 5/7] appease codacy, unused var --- src/pysolid/grid.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/pysolid/grid.py b/src/pysolid/grid.py index 89ceffe..e71f50d 100644 --- a/src/pysolid/grid.py +++ b/src/pysolid/grid.py @@ -11,7 +11,6 @@ # pysolid.calc_solid_earth_tides_grid() -import datetime as dt import os import numpy as np From 630c30ae75572f15808c27db836df73640e28819 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Mon, 17 Apr 2023 11:05:35 -0400 Subject: [PATCH 6/7] return a tuple of arrays from fortran for e/n/u arrays --- src/pysolid/grid.py | 12 ++++-------- src/pysolid/point.py | 9 +++------ src/pysolid/solid.for | 43 +++++++++++++++++++++++++++---------------- 3 files changed, 34 insertions(+), 30 deletions(-) diff --git a/src/pysolid/grid.py b/src/pysolid/grid.py index e71f50d..b81133d 100644 --- a/src/pysolid/grid.py +++ b/src/pysolid/grid.py @@ -73,14 +73,10 @@ def calc_solid_earth_tides_grid(dt_obj, atr, step_size=1e3, display=False, verbo s=(length, width), la=lat_step, lo=lon_step)) ## calc solid Earth tides - fc = solid_grid(dt_obj.year, dt_obj.month, dt_obj.day, - dt_obj.hour, dt_obj.minute, dt_obj.second, - lat0, lat_step, length, - lon0, lon_step, width) - - tide_e = fc[:, :, 0].reshape(length, width) - tide_n = fc[:, :, 1].reshape(length, width) - tide_u = fc[:, :, 2].reshape(length, width) + tide_e, tide_n, tide_u = solid_grid(dt_obj.year, dt_obj.month, dt_obj.day, + dt_obj.hour, dt_obj.minute, dt_obj.second, + lat0, lat_step, length, + lon0, lon_step, width) # resample to the input size if num_step > 1: diff --git a/src/pysolid/point.py b/src/pysolid/point.py index af87fe6..ef1bde3 100644 --- a/src/pysolid/point.py +++ b/src/pysolid/point.py @@ -171,13 +171,10 @@ def calc_solid_earth_tides_point_per_day(lat, lon, date_str, step_sec=60): # calc solid Earth tides t = dt.datetime.strptime(date_str, '%Y%m%d') - fc = solid_point(lat, lon, t.year, t.month, t.day, step_sec) + secs, tide_e, tide_n, tide_u = solid_point( + lat, lon, t.year, t.month, t.day, step_sec + ) - tide_e = fc[:, 1].flatten() - tide_n = fc[:, 2].flatten() - tide_u = fc[:, 3].flatten() - - secs = fc[:, 0].flatten() dt_out = [t + dt.timedelta(seconds=sec) for sec in secs] dt_out = np.array(dt_out) diff --git a/src/pysolid/solid.for b/src/pysolid/solid.for index 7188ccc..06b0a73 100644 --- a/src/pysolid/solid.for +++ b/src/pysolid/solid.for @@ -9,14 +9,15 @@ *** wget -r -l1 --no-parent -R "index.html*" -nH --cut-dirs=3 https://iers-conventions.obspm.fr/content/chapter7/software/dehanttideinel *** Z. Yunjun and S. Sangha, Sep 2020: modify solid() to solid_point/grid() as subroutines. subroutine solid_grid(iyr,imo,idy,ihh,imm,iss, - * glad0,steplat,nlat,glod0,steplon,nlon,output) + * glad0,steplat,nlat,glod0,steplon,nlon,set_east,set_north,set_up) *** calculate solid earth tides (SET) for one spatial grid given the date/time *** Arguments: iyr/imo/idy/ihh/imm/iss - int, date/time for YYYY/MM/DD/HH/MM/SS *** glad0/glad1/steplat - float, north(Y_FIRST)/south/step(negative) in deg *** glod0/glod1/steplon - float, west(X_FIRST) /east /step(positive) in deg -*** Returns: output: 3D array with shape (nlat, nlon, 5) -*** each depth slice: [SET_east, SET_north, SET_up] +*** Returns: set_east - float, east component of SET in m +*** set_north - float, north component of SET in m +*** set_up - float, up component of SET in m implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -24,14 +25,15 @@ integer nlat,nlon double precision glad0,steplat double precision glod0,steplon - real(8), intent(out), dimension(nlat,nlon,3) :: output + real(8), intent(out), dimension(nlat,nlon) :: set_east + real(8), intent(out), dimension(nlat,nlon) :: set_north + real(8), intent(out), dimension(nlat,nlon) :: set_up !***^ leap second table limit flag logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 ! f2py intent(in) iyr,imo,idy,ihh,imm,iss,glad0,steplat,nlat,glod0,steplon,nlon -! f2py dimension(nlat,nlon,3) output -! f2py intent(out) output +! f2py intent(out) set_east,set_north,set_up *** constants @@ -86,8 +88,6 @@ return endif -*** output header - glad1=glad0+nlat*steplat glod1=glod0+nlon*steplon @@ -144,8 +144,10 @@ call mjdciv(mjd,fmjd +0.001d0/86400.d0, * iyr,imo,idy,ihr,imn,sec-0.001d0) - !*** write output to file - output(ilat, ilon, :) = [vt, ut, wt] + !*** write output respective arrays + set_east(ilat, ilon) = vt + set_north(ilat, ilon) = ut + set_up(ilat, ilon) = wt enddo enddo @@ -162,14 +164,17 @@ end *----------------------------------------------------------------------- - subroutine solid_point(glad,glod,iyr,imo,idy,step_sec,output) + subroutine solid_point(glad,glod,iyr,imo,idy,step_sec, + * out_seconds,set_east,set_north,set_up) *** calculate SET at given location for one day with step_sec seconds resolution *** Arguments: glad/glod - float, latitude/longitude in deg *** iyr/imo/idy - int, start date/time in UTC *** step_sec - int, time step in seconds -*** Returns: output -*** each row: [seconds, SET_east, SET_north, SET_up] +*** Returns: out_seconds - float, seconds since start +*** set_east - float, east component of SET +*** set_north - float, north component of SET +*** set_up - float, up component of SET implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -177,13 +182,16 @@ integer iyr,imo,idy integer nloop, step_sec double precision tdel2 - real(8), intent(out), dimension(60*60*24/step_sec,4) :: output + real(8), intent(out), dimension(60*60*24/step_sec) :: out_seconds + real(8), intent(out), dimension(60*60*24/step_sec) :: set_east + real(8), intent(out), dimension(60*60*24/step_sec) :: set_north + real(8), intent(out), dimension(60*60*24/step_sec) :: set_up !*** leap second table limit flag logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 ! f2py intent(in) glad,glod,iyr,imo,idy,step_sec -! f2py intent(out) output +! f2py intent(out) out_seconds,set_east,set_north,set_up *** constants @@ -274,7 +282,10 @@ tsec=ihr*3600.d0+imn*60.d0+sec - output(iloop, :) = [tsec, vt, ut, wt] + out_seconds(iloop) = tsec + set_east(iloop) = vt + set_north(iloop) = ut + set_up(iloop) = wt !*** update fmjd for the next round fmjd=fmjd+tdel2 From ee9705d1f7aa62e2dce54002d5dc36c3dc480f15 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Tue, 18 Apr 2023 12:55:58 +0800 Subject: [PATCH 7/7] solid.for: consistent var name + update comments --- src/pysolid/solid.for | 62 +++++++++++++++++++++---------------------- 1 file changed, 30 insertions(+), 32 deletions(-) diff --git a/src/pysolid/solid.for b/src/pysolid/solid.for index 06b0a73..4f750a3 100644 --- a/src/pysolid/solid.for +++ b/src/pysolid/solid.for @@ -7,17 +7,17 @@ *** J. Gipson and C. Bruyninx. The latest version of dehanttideinel.f and its *** dependent subroutines can be download from IERS conventions website as: *** wget -r -l1 --no-parent -R "index.html*" -nH --cut-dirs=3 https://iers-conventions.obspm.fr/content/chapter7/software/dehanttideinel -*** Z. Yunjun and S. Sangha, Sep 2020: modify solid() to solid_point/grid() as subroutines. +*** Sep 2020: modify solid() to solid_point/grid() as subroutines, Z. Yunjun and S. Sangha. +*** Apr 2023: return numpy arrays instead of writing txt file, S. Staniewicz. + subroutine solid_grid(iyr,imo,idy,ihh,imm,iss, - * glad0,steplat,nlat,glod0,steplon,nlon,set_east,set_north,set_up) + * glad0,steplat,nlat,glod0,steplon,nlon,tide_e,tide_n,tide_u) *** calculate solid earth tides (SET) for one spatial grid given the date/time *** Arguments: iyr/imo/idy/ihh/imm/iss - int, date/time for YYYY/MM/DD/HH/MM/SS *** glad0/glad1/steplat - float, north(Y_FIRST)/south/step(negative) in deg -*** glod0/glod1/steplon - float, west(X_FIRST) /east /step(positive) in deg -*** Returns: set_east - float, east component of SET in m -*** set_north - float, north component of SET in m -*** set_up - float, up component of SET in m +*** glod0/glod1/steplon - float, west (X_FIRST)/east /step(positive) in deg +*** Returns: tide_e/tide_n/tide_u - 2D np.ndarray, east/north/up component of SET in m implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -25,15 +25,15 @@ integer nlat,nlon double precision glad0,steplat double precision glod0,steplon - real(8), intent(out), dimension(nlat,nlon) :: set_east - real(8), intent(out), dimension(nlat,nlon) :: set_north - real(8), intent(out), dimension(nlat,nlon) :: set_up + real(8), intent(out), dimension(nlat,nlon) :: tide_e + real(8), intent(out), dimension(nlat,nlon) :: tide_n + real(8), intent(out), dimension(nlat,nlon) :: tide_u !***^ leap second table limit flag logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 -! f2py intent(in) iyr,imo,idy,ihh,imm,iss,glad0,steplat,nlat,glod0,steplon,nlon -! f2py intent(out) set_east,set_north,set_up + !f2py intent(in) iyr,imo,idy,ihh,imm,iss,glad0,steplat,nlat,glod0,steplon,nlon + !f2py intent(out) tide_e,tide_n,tide_u *** constants @@ -145,9 +145,9 @@ * iyr,imo,idy,ihr,imn,sec-0.001d0) !*** write output respective arrays - set_east(ilat, ilon) = vt - set_north(ilat, ilon) = ut - set_up(ilat, ilon) = wt + tide_e(ilat, ilon) = vt + tide_n(ilat, ilon) = ut + tide_u(ilat, ilon) = wt enddo enddo @@ -165,16 +165,14 @@ *----------------------------------------------------------------------- subroutine solid_point(glad,glod,iyr,imo,idy,step_sec, - * out_seconds,set_east,set_north,set_up) + * secs,tide_e,tide_n,tide_u) *** calculate SET at given location for one day with step_sec seconds resolution -*** Arguments: glad/glod - float, latitude/longitude in deg -*** iyr/imo/idy - int, start date/time in UTC -*** step_sec - int, time step in seconds -*** Returns: out_seconds - float, seconds since start -*** set_east - float, east component of SET -*** set_north - float, north component of SET -*** set_up - float, up component of SET +*** Arguments: glad/glod - float, latitude/longitude in deg +*** iyr/imo/idy - int, start date/time in UTC +*** step_sec - int, time step in seconds +*** Returns: secs - 1D np.ndarray, seconds since start +*** tide_e/tide_n/tide_u - 1D np.ndarray, east/north/up component of SET in m implicit double precision(a-h,o-z) dimension rsun(3),rmoon(3),etide(3),xsta(3) @@ -182,16 +180,16 @@ integer iyr,imo,idy integer nloop, step_sec double precision tdel2 - real(8), intent(out), dimension(60*60*24/step_sec) :: out_seconds - real(8), intent(out), dimension(60*60*24/step_sec) :: set_east - real(8), intent(out), dimension(60*60*24/step_sec) :: set_north - real(8), intent(out), dimension(60*60*24/step_sec) :: set_up + real(8), intent(out), dimension(60*60*24/step_sec) :: secs + real(8), intent(out), dimension(60*60*24/step_sec) :: tide_e + real(8), intent(out), dimension(60*60*24/step_sec) :: tide_n + real(8), intent(out), dimension(60*60*24/step_sec) :: tide_u !*** leap second table limit flag logical lflag common/stuff/rad,pi,pi2 common/comgrs/a,e2 -! f2py intent(in) glad,glod,iyr,imo,idy,step_sec -! f2py intent(out) out_seconds,set_east,set_north,set_up + !f2py intent(in) glad,glod,iyr,imo,idy,step_sec + !f2py intent(out) secs,tide_e,tide_n,tide_u *** constants @@ -282,10 +280,10 @@ tsec=ihr*3600.d0+imn*60.d0+sec - out_seconds(iloop) = tsec - set_east(iloop) = vt - set_north(iloop) = ut - set_up(iloop) = wt + secs(iloop) = tsec + tide_e(iloop) = vt + tide_n(iloop) = ut + tide_u(iloop) = wt !*** update fmjd for the next round fmjd=fmjd+tdel2