Skip to content

Commit

Permalink
update solid.for to return arrays instead of using text file (#56)
Browse files Browse the repository at this point in the history 
+ solid.for/solid_grid/point(): return numpy arrays instead of writing text file, by leveraging the `!f2py intent(in/out)` comment syntax. This brings two improvements:
   1. greatly improve the communication between the fortran and python code, thus, the flexibility of pysolid
   2. >2X speedup of point SET calculation

+ update `solid_grid/point()` usage in grid/point.py

---------

Co-authored-by: Zhang Yunjun <[email protected]>
  • Loading branch information
@scottstanie @yunjunz
scottstanie and yunjunz authored Apr 18, 2023
1 parent 94bec28 commit 18dcc0f
Show file tree
Hide file tree
Showing 4 changed files with 84 additions and 142 deletions.
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +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"])
],
)
31 changes: 5 additions & 26 deletions src/pysolid/grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@


import os
import tempfile

import numpy as np
from skimage.transform import resize
Expand Down Expand Up @@ -73,31 +72,11 @@ 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
with tempfile.NamedTemporaryFile(prefix="pysolid_", suffix=".txt") as fp:
vprint('SOLID : calculating / writing data to txt file: {}'.format(fp.name))

# Run twice to circumvent fortran bug which cuts off last file in loop
# - Simran, Jun 2020
for _ in range(2):
solid_grid(fp.name, 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(fp.name))
grid_size = int(length * width)
fc = np.loadtxt(fp.name,
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)
## calc solid Earth tides
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:
Expand Down
30 changes: 6 additions & 24 deletions src/pysolid/point.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -170,29 +169,12 @@ 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

# create a temporary text file so it doesn't get overwritten by competing processes
with tempfile.NamedTemporaryFile(prefix="pysolid_", suffix=".txt") as fp:
# 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(fp.name, 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(fp.name,
dtype=float,
delimiter=',',
skiprows=0,
max_rows=num_row)

tide_e = fc[:, 1].flatten()
tide_n = fc[:, 2].flatten()
tide_u = fc[:, 3].flatten()

secs = fc[:, 0].flatten()
# calc solid Earth tides
t = dt.datetime.strptime(date_str, '%Y%m%d')
secs, tide_e, tide_n, tide_u = solid_point(
lat, lon, t.year, t.month, t.day, step_sec
)

dt_out = [t + dt.timedelta(seconds=sec) for sec in secs]
dt_out = np.array(dt_out)

Expand Down
163 changes: 72 additions & 91 deletions src/pysolid/solid.for
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,36 +7,33 @@
*** 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.

subroutine solid_grid(txt_file,iyr,imo,idy,ihh,imm,iss,
* glad0,steplat,nlat,
* glod0,steplon,nlon)
*** 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,tide_e,tide_n,tide_u)

*** 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: latitude, longitude, SET_east, SET_north, SET_up
*** 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)
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
double precision glad0,steplat
double precision glod0,steplon
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

*** open output file

lout=1
open(lout,file=txt_file,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 intent(out) tide_e,tide_n,tide_u

*** constants

Expand All @@ -52,63 +49,55 @@
*** 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

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)

Expand Down Expand Up @@ -155,53 +144,52 @@
call mjdciv(mjd,fmjd +0.001d0/86400.d0,
* iyr,imo,idy,ihr,imn,sec-0.001d0)

!*** write output to file
write(lout,'(*(G0.9,:,", "))') glad,glod,vt,ut,wt
!*** write output respective arrays
tide_e(ilat, ilon) = vt
tide_n(ilat, ilon) = ut
tide_u(ilat, ilon) = wt

enddo
enddo

*** 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(txt_file,glad,glod,iyr,imo,idy,step_sec)
subroutine solid_point(glad,glod,iyr,imo,idy,step_sec,
* secs,tide_e,tide_n,tide_u)

*** 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
*** iyr/imo/idy - int, start date/time in UTC
*** step_sec - int, time step in seconds
*** Returns: seconds, SET_east, SET_north, SET_up
*** 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)
character(len=*), intent(in) :: txt_file
dimension rsun(3),rmoon(3),etide(3),xsta(3)
double precision glad,glod
integer iyr,imo,idy
integer nloop, step_sec
double precision tdel2
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

*** open output file

lout=1
open(lout,file=txt_file,form='formatted',status='unknown')
write(lout,'(a)') '# program solid -- UTC version -- 2018jun01'
!f2py intent(in) glad,glod,iyr,imo,idy,step_sec
!f2py intent(out) secs,tide_e,tide_n,tide_u

*** constants

Expand All @@ -217,34 +205,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)

Expand Down Expand Up @@ -272,11 +256,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
Expand All @@ -296,9 +278,12 @@

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

secs(iloop) = tsec
tide_e(iloop) = vt
tide_n(iloop) = ut
tide_u(iloop) = wt

!*** update fmjd for the next round
fmjd=fmjd+tdel2
Expand All @@ -309,16 +294,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)
Expand Down

0 comments on commit 18dcc0f

Please sign in to comment.