pyacs.lib package¶

Subpackages¶

pyacs.lib.faultslip package

Submodules¶

pyacs.lib.astrotime module¶

AstroTime.py is a traduction from perl module Astro::Time distributed by CPAN and written by Chris Phillips (Chris.Phillips@csiro.au). Conversion to python made by Jean-Mathieu Nocquet (Geoazur - CNRS-IRD-OCA-Univ. Nice, nocquet@geoazur.unice.fr)

AstroTime contains a set of Python routines for time based conversions, such as conversion between calendar dates and Modified Julian day and conversion from UT to local sidereal time. Included are routines for conversion between numerical and string representation of angles. All string functions removed.

Conversion from and to datetime objects have also been added for convenience.

note: This module is intended dates manipulation, but leap seconds for instance are not handled. See the GPSTime module for GPS, UTC times conversions.
note: python 3.6 : all / operator changed to // to force integer operations
warning: !!! all conversions not specifying ut or uts use 12h00mn00s as a default. ut is the decimal fraction of day and uts is the number of seconds since 00h00m00.0s

pyacs.lib.astrotime.cal2dayno(day, month, year)[source]¶

Returns the day of year (doy).

Parameters: day,month,year –
Type: int,int,int
Returns: doy of year as int
Example

>>> import pyacs
>>> pyacs.cal2dayno(1,1,2000)
>>> 1
>>> pyacs.cal2dayno(29,2,2000)
>>> 60
>>> pyacs.cal2dayno(29.,2,2000)
>>> 60.
>>> pyacs.cal2dayno(29,2,2001)
>>> !!!  29  out of range
>>> !!! bad day

pyacs.lib.astrotime.cal2mjd(day, month, year, ut=0.5)[source]¶

Converts a calendar date (universal time) into modified Julian day number. mjd Modified Julian day (JD-2400000.5)

Parameters: month, year (day,) –

:param ut : day fraction ([0.,1.[) :returns: mjd (modified Julian day (julian day -2400000.5)) :rtype: float :note: Be aware that when ut is not provided, then the middle of the day is used. See example.

Example

>>> import pyacs
>>> pyacs.cal2mjd(29,2,2000)
>>> 51603.5
>>> pyacs.cal2mjd(29,2,2000,ut=0.0)
>>> 51603.0

pyacs.lib.astrotime.cal2decyear(day, month, year, ut=0.5)[source]¶

converts a calendar date to decimal year

Parameters: month, year (day,) –

:param ut : day fraction ([0.,1.[) :returns: decimal year :rtype: float :note: Be aware that when ut is not provided, then the middle of the day is used.

pyacs.lib.astrotime.cal2datetime(day, month, year, ut=0.5)[source]¶

converts a calendar date to a datime python instance

Parameters: month, year (day,) –

:param ut : day fraction ([0.,1.[), optional

:note:!!! the returned datetime is at 12:00, that is ut=0.5 by default

pyacs.lib.astrotime.dayno2cal(dayno, year)[source]¶

Returns the day and month corresponding to dayno of year.

Parameters: year (dayno,) –
Returns: day, month

pyacs.lib.astrotime.dayno2mjd(dayno, year, ut=0.5)[source]¶

converts a dayno and year to modified Julian day

Parameters

year (dayno,) –
ut – day fraction, optional

Returns

modified Julian day (julian day - 2400000.5)

Note

Be aware that when ut is not provided, then the middle of the day is used. See example.

Example

>>> import pyacs
>>> pyacs.dayno2mjd(60,2000)
>>> 51603.5
>>> pyacs.dayno2mjd(60.,2000,ut=0)
>>> 51603.0

pyacs.lib.astrotime.dayno2decyear(dayno, year, ut=0.5)[source]¶

converts julian day to decimal year

Parameters

year (dayno,) –
ut – day fraction, optional

Returns

decimal year

Note

!!! accounts for leap years with 366 days

Note

!!! unless ut is specified, the returned decimal years is at 12:00, that is ut=0.5 by default

Example

>>> import pyacs
>>> pyacs.dayno2decyear(60,1999)
>>> 1999.16301369863
>>> pyacs.dayno2decyear(60,2000)
>>> 2000.1625683060108
>>> pyacs.dayno2decyear(60,2001)
>>> 2001.16301369863
>>> pyacs.dayno2decyear(60,2001,ut=0.)
>>> 2001.1616438356164

pyacs.lib.astrotime.dayno2datetime(dayno, year, ut=0.5)[source]¶

julian day to datetime python instance

Parameters

year (dayno,) –
ut – day fraction, optional

Returns

datetime

:note:!!! the returned datetime is at 12:00, that is uts=24.0*60.0*60/2.0 by default

:note:!!! specify uts for other time of the day

pyacs.lib.astrotime.mjd2cal(mjd)[source]¶

Converts a modified Julian day number into calendar date (universal time). (based on the slalib routine sla_djcl).

Parameters: mjd – modified Julian day (JD-2400000.5)
Returns: day,month,year,ut
Return type: int,int,int,float
Note: ut is the the day fraction in [0., 1.[.

pyacs.lib.astrotime.mjd2dayno(mjd)[source]¶

converts a modified Julian day into year and dayno (universal time).

Parameters: mjd – modified julian day
Returns: dayno,year,ut
Return type: int,int,float
Example

>>> import pyacs
>>> pyacs.mjd2dayno(51603.0)
>>> (60, 2000, 0.0)
>>> pyacs.mjd2dayno(51603.5)
>>> (60, 2000, 0.5)

pyacs.lib.astrotime.mjd2decyear(mjd)[source]¶

converts a modified julian day to decimal year conversion

Parameters: mjd – modified julian day
Returns: decimal year

pyacs.lib.astrotime.mjd2datetime(mjd)[source]¶

converts Modified Julian Day date to datime python instance

Parameters: mjd – modified julian day
Returns: datetime instance

pyacs.lib.astrotime.mjd2gpsweek(mjd)[source]¶

converts a modified julian day to gps week and gps day of week

Parameters: mjd – modified julian day
Returns: gps_week,day_of_week

pyacs.lib.astrotime.decyear2cal(decyear)[source]¶

decimal year to calendar date conversion

Parameters: decyear – decimal year
Returns: day of month,month,ut
Note: the input decimal year is assumed to account for leap years, that is the day of year is the decimal of decyear * number of days in the year.

pyacs.lib.astrotime.decyear2dayno(decyear)[source]¶

converts decimal year to the day of year

Parameters: decyear – decimal year
Returns: day of year,ut
Note: the input decimal year is assumed to account for leap years, that is the day of year is the decimal of decyear * number of days in the year.

pyacs.lib.astrotime.decyear2mjd(decyear)[source]¶

decimal year to modified julian day conversion

Parameters: decyear – decimal year
Returns: modified julian day

pyacs.lib.astrotime.decyear2datetime(decyear)[source]¶

converts a decimal date to a datime python instance

Parameters: decyear – decimal year
Returns: datetime instance

pyacs.lib.astrotime.decyear2epoch(decyear)[source]¶

converts decimal year to formatted epoch %02d:%03d:%05d

Parameters: decyear – decimal year

:returns : ‘yy:doy:sec’

Note: the input decimal year is assumed to account for leap years, that is the day of year is the decimal of decyear * number of days in the year.
Example

>>> import pyacs
>>> pyacs.cal2decyear(1,3,2001)
>>> 2001.16301369863
>>> pyacs.dec
>>> pyacs.decyear2epoch(2001.16301369863)
>>> '01:060:43200'
>>> pyacs.decyear2epoch(2000.16301369863)
>>> '00:060:57284'
>>> pyacs.decyear2epoch(pyacs.cal2decyear(1,3,2000))
>>> '00:061:43200'

pyacs.lib.astrotime.datetime_from_calarray(calarray, hour=12, minute=0, sec=0, ms=0)[source]¶

Build a datetime array from an array of calendar dates

Parameters: calarray – 1D or 2D numpy array. Columns order are year,month,mday,hour,minute,second,microsecond
Returns: numpy 1D array of datime instances
Note: year,month,mday are mandatory. If not provided then hour=12, minute=0, sec=0,ms=0 (middle of the day)

pyacs.lib.astrotime.datetime2cal(datetime)[source]¶

converts from python datetime to calendar date

Parameters: datetime – datetime instance
Returns: year,month,monthday,ut

pyacs.lib.astrotime.datetime2dayno(datetime)[source]¶

Converts python datetime to dayno,year,ut

Parameters: datetime – pythob datetime instance
Returns: year,dayno,ut

pyacs.lib.astrotime.datetime2mjd(datetime)[source]¶

converts a python datetime instance to a modified julian day

Parameters: datetime – python datetime instance
Returns: modified julian day

pyacs.lib.astrotime.datetime2decyear(datetime)[source]¶

converts a python datetime instance to decimal year

Parameters: datetime – python datetime instance or array of it
Returns: decimal year

pyacs.lib.astrotime.datetime_round_second(my_datetime)[source]¶

rounds a python datetime instance to the nearest second

Parameters: datetime – python datetime instance or array of it
Returns: new rounded datetime

pyacs.lib.astrotime.gpsweek2mjd(gpsweek, dow)[source]¶

converts a GPS week and dow to modified julian day

Parameters: gpsweek,dow – gps week, day of week (0 is sunday, 6 saturday)
Returns: modified julian day

pyacs.lib.astrotime.epoch2decyear(epoch)[source]¶

Converts an epoch string used in the SINEX format ‘YY:DOY:SOD’ e.g.’17:100:43185’ into a decimal year

Parameters: epoch – epoch formatted string or list of epochs
Returns: decimal year or array of decimal years

pyacs.lib.astrotime.jd2mjd(jd)[source]¶

Converts a Julian day to a Modified Julian day

Parameters: jd – Julian day
Returns: jd-2400000.5

pyacs.lib.astrotime.mjd2jd(mjd)[source]¶

converts a Modified Julian day to Julian day

Parameters: mjd – Modified Julian day
Returns: mjd + 2400000.5

pyacs.lib.astrotime.yr2year(yr)[source]¶

converts two digits year to four digits year

Parameters: yr – 2-digits year (e.g. 01)
Returns: 4-digit years (2001). This is yr+2000 if yr<80 yr+1900 otherwise.
Note: This is an heritage from a bad practice in the geodesy community assuming that origin of the universe started with GPS. See example.
Example

>>> import pyacs
>>> pyacs.yr2year(96)
>>> 1996
>>> pyacs.yr2year(11)
>>> 2011

pyacs.lib.astrotime.year2yr(year)[source]¶

converts 4-digits year to 2-digits year

Parameters: yr – 4-digits year (e.g. 2001)
Returns: 2-digit years (01). This is yr-2000 if yr<80 yr-1900 otherwise.
Note: This is an heritage from a bad practice in the geodesy community assuming that origin of the universe started with GPS in 1980. See example.
Example

>>> import pyacs
>>> pyacs.year2yr(1996)
>>> 96
>>> pyacs.year2yr(2011)
>>> 11

pyacs.lib.astrotime.leap_year(year)[source]¶

Returns true if year is a leap year.

Parameters: year – year in YYYY
Returns: True if year is leap, False otherwise

pyacs.lib.astrotime.uts2hmsmicros(uts)[source]¶

converts decimal seconds on a day (0-86400) to hours, min, seconde, microseconds

Parameters: uts – decimal seconds since 00:00:00.0. Must be in the [0.,86400.[ range.
Returns: hours,minutes,seconds,microseconds
Return type: int,int,int,float
Example

>>> import pyacs
>>> pyacs.uts2hmsmicros(86399.999999)
>>> (23, 59, 59, 999999.0000069374)
>>> pyacs.uts2hmsmicros(86400)
>>> !!!  86400 uts out of range [0- 86400.0 [

pyacs.lib.astrotime.hmsmicros2uts(h, mn, s=0, microsecond=0)[source]¶

converts hour minute seconds to uts

:param hours,minutes,seconds,microseconds :returns: uts

pyacs.lib.astrotime.hmsmicros2ut(h, mn, s=0, microsecond=0)[source]¶

converts hour minute seconds to fractional day

:param hours,minutes,seconds,microseconds :returns: ut

pyacs.lib.astrotime.ut2uts(ut)[source]¶

converts fractional day to seconds

Parameters: ut – fractional day in [0.,1.[
Returns: ut * 60.*60.*24.

pyacs.lib.astrotime.uts2ut(uts)[source]¶

converts uts in seconds to fractional day

Parameters: uts – seconds in [0.,86400.[
Returns: uts / (60.*60.*24.)

pyacs.lib.astrotime.day_since_decyear(decyear, ref_decyear)[source]¶

Calculates the number of days since a reference date.

Parameters: decyear – decimal year
Returns: days elapsed since decyear (float)
Note: negative number means before the reference date
Note: a useful function for postseismic analysis.
Example

>>> import pyacs
>>> ref_decyear=pyacs.cal2decyear(16,4,2016,ut=pyacs.hmsmicros2ut(23, 58, 33)) # Ecuador, Mw 7.8 EQ
>>> pyacs.day_since_decyear(pyacs.cal2decyear(17,4,2016,ut=pyacs.hmsmicros2ut(23, 58, 32)),ref_decyear)
>>> 0.9999884259232203
>>> pyacs.day_since_decyear(pyacs.cal2decyear(17,4,2016,ut=pyacs.hmsmicros2ut(23, 58, 33)),ref_decyear)
>>> 0.9999999999854481
>>> pyacs.day_since_decyear(pyacs.cal2decyear(16,4,2017,ut=pyacs.hmsmicros2ut(23, 58, 33)),ref_decyear)
>>> 364.9999999999345
>>> ld=pyacs.mjd2decyear(np.arange(pyacs.cal2mjd(13,4,2016),pyacs.cal2mjd(19,4,2016)))
>>> ld
>>> array([ 2016.283,  2016.286,  2016.288,  2016.291,  2016.294,  2016.296])
>>> pyacs.day_since_decyear(ld,ref_decyear)
>>> array([-3.499, -2.499, -1.499, -0.499,  0.501,  1.501])

pyacs.lib.astrotime.guess_date(date)[source]¶

Tries to guess a date from input parameter returns a decimal year

Parameters: date – date as string or float
Returns: decimal year
Example

>>> import pyacs
>>> pyacs.guess_date('2016/04/16')
>>> 2016.2909836065573

pyacs.lib.astrotime.decyear2seconds(np_decyear, rounding='day')[source]¶

converts decimal year to a seconds since 1980/1/1/0/0/0

Parameters: rounding – controls rounding. ‘day’ means at 12:00:00,

‘hour’ at 00:00, ‘minute’ at the current minute with 00 seconds, ‘second’ at the integer of the current second. :return: 1D integer (np.int64) numpy array

pyacs.lib.astrotime.datetime2seconds(np_datetime)[source]¶: converts dattime to seconds since 1980/1/1/0/0/0 :param datetime: :return: 1D numpy array of seconds

pyacs.lib.astrotime.seconds2datetime(seconds)[source]¶

converts seconds since 1980/1/1/0/0/0 to a datetime object

:return:a 1D numpy array of datetime object

pyacs.lib.coordinates module¶

Coordinates library.

– Local/Geocentric frame conversion

– Geodetic/Geocentric frame conversion

– Spherical/Geocentric conversion

– Geodetic/Flat Earth conversion

pyacs.lib.coordinates.mat_rot_general_to_local(lam, phi, unit='radians')[source]¶

Generates the conversion matrix R from general geocentric cartesian coordinates (XYZ) to local cartesian coordinates (ENU):

Parameters

lam,phi – longitude
unit – ‘radians’ or ‘dec_deg’ (default is ‘radians’)

Returns

R as a 2D numpy array

Note

R = [[-sin(lambda) cos(lambda) 0 ],

[-sin(phi)*cos(lambda) , -sin(phi)*sin(lamda) , cos(phi)],

[ cos(phi)*cos(lambda) , cos(phi)*sin(lamda) , sin(phi)]]

pyacs.lib.coordinates.mat_rot_local_to_general(lam, phi, unit='radians')[source]¶

Generates the conversion matrix R from local cartesian coordinates (ENU) to general geocentric cartesian coordinates (XYZ):

Parameters

lam,phi – longitude in radians
unit – ‘radians’ or ‘dec_deg’ (default is ‘radians’)

Returns

R as a 2D numpy array

Note

Since R is orthogonal, it is the inverse and also the transpose of the conversion matrix from general geocentric cartesian coordinates (XYZ) to local cartesian coordinates (ENU)

pyacs.lib.coordinates.denu_at_x0y0z0_to_xyz(de, dn, du, x0, y0, z0)[source]¶

Converts a local vector [dn,de,du] at [x0,y0,z0] in geocentric cartesian coordinates into X,Y,Z geocentric cartesian coordinates

Parameters

de,dn,du – east,north, up components of the local vector in meters
x0,x0,z0 – reference point for the local vector in geocentric cartesian coordinates

Returns

x,y,z in geocentric cartesian coordinates

pyacs.lib.coordinates.xyz2geo(x, y, z, A=6378137.0, E2=0.006694380022903, unit='radians')[source]¶

Converts geocentric cartesian coordinates (XYZ) to geodetic coordinates (lambda,phi,h)

Parameters

x,y,z – XYZ coordinates in meters
unit – ‘radians’ or ‘dec_deg’ (default is ‘radians’)

Returns

long,lat,he: longitude and latitude in radians, he in meters above the ellipsoid

Note

Default ellipsoid is GRS80 used for WGS84 with

A = 6378137. # semi major axis = equatorial radius

E2 = 0.00669438003 # eccentricity and then

F = 1.0 - sqrt(1-E2) # flattening

Ref: Bowring, 1985, The accuracy of geodetic latitude and height equations, Survey Review, 28, 202-206.

pyacs.lib.coordinates.geo2xyz(llambda, phi, he, unit='radians', A=6378137.0, E2=0.006694380022903)[source]¶

Converts geodetic coordinates (long,lat,he) to geocentric cartesian coordinates (XYZ)

Parameters

llambda,phi – longitude, latitude
he – ellipsoidal height in meters
unit – ‘radians’ or ‘dec_deg’ for llambda and phi (default is ‘radians’)

Returns

x,y,z in meters

Note

Default ellipsoid is GRS80 used for WGS84 with

A = 6378137. # semi major axis = equatorial radius

E2 = 0.00669438003 # eccentricity and then

F = 1.0 - sqrt(1-E2) # flattening

pyacs.lib.coordinates.wnorm(phi, A=6378137.0, E2=0.006694380022903)[source]¶: Calculates the geodetic radius normal to the ellipsoid

pyacs.lib.coordinates.xyz2geospheric(x, y, z, unit='radians')[source]¶

Converts geocentric cartesian coordinates (XYZ) to geo-spherical coordinates (longitude,latitude,radius).

Parameters

x,y,z – geocentric cartesian coordinates in meters
unit – ‘radians’ or ‘dec_deg’ (default is ‘radians’)

Returns

longitude,latitude (in radians), radius from the Earth’s center in meters

Note

Be aware that the obtained coordinates are not what is usually taken as spherical coordinates, which uses co-latitude

pyacs.lib.coordinates.xyz_spherical_distance(x1, y1, z1, x2, y2, z2, Rt=6371000.0)[source]¶

Returns the spherical distance between two points of XYZ coordinates x1,y1,z1 and x2,y2,z2

Parameters

x1,y1,z1,x2,y2,z2 – in meters
Rt – mean Earth’s radius in meters (default 6 371 000.0 m)

Returns

distance in meters

pyacs.lib.coordinates.geo_spherical_distance(lon1, lat1, h1, lon2, lat2, h2, unit='radians', Rt=6371000.0)[source]¶

Returns the spherical distance between two points of geodetic coordinates lon1,lat1,lon2,lat2

Parameters

lon1,lat1,h1,lon2,lat2,h2 – longitude, latitude, radius
unit – ‘radians’ or ‘dec_deg’ for lon & lat (default is ‘radians’)
Rt – mean Earth’s radius in meters (default 6 371 000.0 m)

Returns

distance in meters

pyacs.lib.coordinates.azimuth_to_en(azimuth)[source]¶

converts an azimuth (clockwise from north) to east and north unit vector

Parameters: azimuth – azimuth in decimal degrees
Return east,north: unit vector
Note: calculation is made using the sphere

pyacs.lib.coordinates.geo2flat_earth(longitude, latitude)[source]¶

Converts geographical coordinates to flat earth approximation uses pyproj and web Mercator projection.

Parameters: longitude,latitude – geographical (ellipsoid) coordinates in decimal degrees. Also works with 1D numpy arrays
Returns x,y: in km

pyacs.lib.coordinates.flat_earth2geo(x, y)[source]¶

Converts web Mercator coordinates (units km) to geographical coordinates uses pyproj and web Mercator projection.

Parameters: longitude,latitude – geographical (ellipsoid) coordinates in decimal degrees. Also works with 1D numpy arrays
Returns x,y: in km

pyacs.lib.coordinates.spherical_baseline_length_rate(slon, slat, sve, svn, elon, elat, eve, evn, sigma=None, verbose=False)[source]¶

calculates the baseline (great circle) length rate of change between two points

Parameters

slon,slat – coordinates of profile start point (decimal degrees)
sve,svn – east and north components of start point (mm/yr)
elon,elat – coordinates of profile end point (decimal degrees)
eve,evn – east and north components of end point (mm/yr)
sigma – list of velocities uncertainties: [sig_sve, sig_svn, corr_sven, sig_eve, sig_evn, corr_even]

:return : length, rate, 1D strain rate

pyacs.lib.errors module¶

Base class for exceptions in PYACS lib core module

exception pyacs.lib.errors.PyacsError[source]¶

Bases: BaseException

Base class for exceptions in PYACS lib core module

exception pyacs.lib.errors.BadOption[source]¶

Bases: Exception

Exception raised for an option not allowed.

exception pyacs.lib.errors.ReadingFile(file_name)[source]¶: Bases: pyacs.lib.errors.PyacsError

pyacs.lib.euclid module¶

euclid graphics maths module

Documentation and tests are included in the file “euclid.rst”, or online at https://github.com/ezag/pyeuclid/blob/master/euclid.rst

class pyacs.lib.euclid.Vector2(x=0, y=0)[source]¶

Bases: object

x¶

y¶

copy()¶

magnitude()¶

magnitude_squared()[source]¶

normalize()[source]¶

normalized()[source]¶

dot(other)[source]¶

cross()[source]¶

reflect(normal)[source]¶

angle(other)[source]¶: Return the angle to the vector other

project(other)[source]¶: Return one vector projected on the vector other

class pyacs.lib.euclid.Vector3(x=0, y=0, z=0)[source]¶

Bases: object

x¶

y¶

z¶

copy()¶

magnitude()¶

magnitude_squared()[source]¶

normalize()[source]¶

normalized()[source]¶

dot(other)[source]¶

cross(other)[source]¶

reflect(normal)[source]¶

rotate_around(axis, theta)[source]¶: Return the vector rotated around axis through angle theta. Right hand rule applies

angle(other)[source]¶: Return the angle to the vector other

project(other)[source]¶: Return one vector projected on the vector other

class pyacs.lib.euclid.Matrix3[source]¶

Bases: object

copy()¶

identity()[source]¶

scale(x, y)[source]¶

translate(x, y)[source]¶

rotate(angle)[source]¶

classmethod new_identity()[source]¶

classmethod new_scale(x, y)[source]¶

classmethod new_translate(x, y)[source]¶

classmethod new_rotate(angle)[source]¶

determinant()[source]¶

inverse()[source]¶

a¶

b¶

c¶

e¶

f¶

g¶

i¶

j¶

k¶

class pyacs.lib.euclid.Matrix4[source]¶

Bases: object

copy()¶

transform(other)[source]¶

identity()[source]¶

scale(x, y, z)[source]¶

translate(x, y, z)[source]¶

rotatex(angle)[source]¶

rotatey(angle)[source]¶

rotatez(angle)[source]¶

rotate_axis(angle, axis)[source]¶

rotate_euler(heading, attitude, bank)[source]¶

rotate_triple_axis(x, y, z)[source]¶

transpose()[source]¶

transposed()[source]¶

classmethod new(*values)[source]¶

classmethod new_identity()[source]¶

classmethod new_scale(x, y, z)[source]¶

classmethod new_translate(x, y, z)[source]¶

classmethod new_rotatex(angle)[source]¶

classmethod new_rotatey(angle)[source]¶

classmethod new_rotatez(angle)[source]¶

classmethod new_rotate_axis(angle, axis)[source]¶

classmethod new_rotate_euler(heading, attitude, bank)[source]¶

classmethod new_rotate_triple_axis(x, y, z)[source]¶

classmethod new_look_at(eye, at, up)[source]¶

classmethod new_perspective(fov_y, aspect, near, far)[source]¶

determinant()[source]¶

inverse()[source]¶

get_quaternion()[source]¶: Returns a quaternion representing the rotation part of the matrix. Taken from: http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q55

a¶

b¶

c¶

d¶

e¶

f¶

g¶

h¶

i¶

j¶

k¶

l¶

m¶

n¶

o¶

p¶

class pyacs.lib.euclid.Quaternion(w=1, x=0, y=0, z=0)[source]¶

Bases: object

w¶

x¶

y¶

z¶

copy()¶

magnitude()¶

magnitude_squared()[source]¶

identity()[source]¶

rotate_axis(angle, axis)[source]¶

rotate_euler(heading, attitude, bank)[source]¶

rotate_matrix(m)[source]¶

conjugated()[source]¶

normalize()[source]¶

normalized()[source]¶

get_angle_axis()[source]¶

get_euler()[source]¶

get_matrix()[source]¶

classmethod new_identity()[source]¶

classmethod new_rotate_axis(angle, axis)[source]¶

classmethod new_rotate_euler(heading, attitude, bank)[source]¶

classmethod new_rotate_matrix(m)[source]¶

classmethod new_interpolate(q1, q2, t)[source]¶

class pyacs.lib.euclid.Geometry[source]¶

Bases: object

intersect(other)[source]¶

connect(other)[source]¶

distance(other)[source]¶

class pyacs.lib.euclid.Point2(x=0, y=0)[source]¶

Bases: pyacs.lib.euclid.Vector2, pyacs.lib.euclid.Geometry

intersect(other)[source]¶

connect(other)[source]¶

x¶

y¶

class pyacs.lib.euclid.Line2(*args)[source]¶

Bases: pyacs.lib.euclid.Geometry

p¶

v¶

copy()¶

property p1¶

property p2¶

intersect(other)[source]¶

connect(other)[source]¶

class pyacs.lib.euclid.Ray2(*args)[source]¶

Bases: pyacs.lib.euclid.Line2

p¶

v¶

class pyacs.lib.euclid.LineSegment2(*args)[source]¶

Bases: pyacs.lib.euclid.Line2

magnitude_squared()[source]¶

property length¶

p¶

v¶

class pyacs.lib.euclid.Circle(center, radius)[source]¶

Bases: pyacs.lib.euclid.Geometry

c¶

r¶

copy()¶

intersect(other)[source]¶

connect(other)[source]¶

tangent_points(p)[source]¶

class pyacs.lib.euclid.Point3(x=0, y=0, z=0)[source]¶

Bases: pyacs.lib.euclid.Vector3, pyacs.lib.euclid.Geometry

intersect(other)[source]¶

connect(other)[source]¶

x¶

y¶

z¶

class pyacs.lib.euclid.Line3(*args)[source]¶

Bases: object

p¶

v¶

copy()¶

property p1¶

property p2¶

intersect(other)[source]¶

connect(other)[source]¶

class pyacs.lib.euclid.Ray3(*args)[source]¶

Bases: pyacs.lib.euclid.Line3

p¶

v¶

class pyacs.lib.euclid.LineSegment3(*args)[source]¶

Bases: pyacs.lib.euclid.Line3

magnitude_squared()[source]¶

property length¶

p¶

v¶

class pyacs.lib.euclid.Sphere(center, radius)[source]¶

Bases: object

c¶

r¶

copy()¶

intersect(other)[source]¶

connect(other)[source]¶

class pyacs.lib.euclid.Plane(*args)[source]¶

Bases: object

n¶

k¶

copy()¶

intersect(other)[source]¶

connect(other)[source]¶

pyacs.lib.euler module¶

Euler poles manipulation

pyacs.lib.euler.rot2euler(Wx, Wy, Wz)[source]¶

Converts a rotation rate vector Wx Wy Wz in radians/yr into an Euler pole (long. lat. deg/Myr) into geographical (spherical) coordinates and angular velocity

Parameters: Wx,Wy,Wz – rotation rate vector in geocentric cartesian coordinates with units of radians/yr
Returns longitude,latitude,omega: longitude and latitude of the Euler pole in radians and the

angular velocity in decimal degrees per Myr.

Note: longitude and latitude are relative to the sphere, not the ellipsoid. This is because

Euler pole and rigid rotations only have sense on a sphere.

pyacs.lib.euler.euler2rot(lon, lat, omega)[source]¶

Converts Euler poles (long., lat., deg/Myr) into cartesian geocentric: rotation rate vector Wx Wy Wz in radians/yr

Parameters longitude,latitude,omega: longitude and latitude of the Euler pole in decimal degrees and the

angular velocity in decimal degrees per Myr. :returns Wx Wy Wz: in radians/yr

Note: longitude and latitude are relative to the sphere, not the ellipsoid.

pyacs.lib.euler.euler_uncertainty(w, vcv)[source]¶

Calculates Euler pole parameters uncertainty

Parameters

w – rotation vector in XYZ coordinates as numpy 1D array
vcv – covariance of w as numpy 2D array

Returns

vcv_euler as numpy 2D array

pyacs.lib.euler.vel_from_euler(lonp, latp, lon_euler, lat_euler, omega_euler)[source]¶

Return the horizontal velocity predicted at lonp,latp from an Euler pole

Parameters

lonp,latp – longitude,latitude in decimal degrees where velocity will be predicted
lon_euler,lat_euler_omega_euler – longitude and latitude of the Euler pole in decimal degrees and the

angular velocity in decimal degrees per Myr.

Returns: ve,vn in mm/yr

pyacs.lib.euler.pole_matrix(coor)[source]¶

Calculates the matrix relating the horizontal velocity to a rotation rate vector. Given a 2D-numpy array of n positions [lonp , latp] in decimal degrees the return matrix is W so that np.dot( W , w ) gives a 2D-numpy array of [ve1,vn1,ve2,vn2,….] expressed in m/yr for a rotation rate vector in rad/yr

Parameters: coor – 2D numpy array of [lon, lat] in decimal degrees
Returns: the pole matrix as a 2D-numpy array

pyacs.lib.euler.pole_matrix_fault(coor, strike, order=None)[source]¶

Calculates the matrix relating the along strike and normal slip components of a fault to a rotation rate vector. Given a 2D-numpy array of n positions [lonp , latp] in decimal degrees and strike counter-clockwise from north the return matrix is W so that np.dot( W , w ) gives a 2D-numpy array of [ss1,ns1,ss2,ns2,….] expressed in m/yr for a rotation rate vector in rad/yr

:param coor : 2D numpy array of [lon, lat] in decimal degrees :param strike: 1D numpy array of [strike] in decimal degrees

Returns: the pole matrix as a 2D-numpy array
Note: np.dot( W , w ).reshape(-1,2) gives the along strike,and normal components in two columns
Note: the value are given for the hanging-wall block (right-polygon)

pyacs.lib.glinalg module¶

Linear algebra for Geodesy problems

pyacs.lib.glinalg.ls(G, d, verbose=False)[source]¶

Solve the least-squares (LS) problem m so that |Gx-d|**2 is minimum.

Parameters

G – m x n model matrix as 2D numpy array
d – m 1D numpy observation vector
verbose – verbose mode

Returns

x,chi2: m (1D numpy array of dim m), chi2 (chi-square)

Note

solved through numpy.linalg.lstsq

pyacs.lib.glinalg.lsw(G, d, std)[source]¶

Solve the least-squares (LS) with data uncertainties provided as a vector

Parameters

G – m x n model matrix as 2D numpy array
d – m 1D numpy observation vector
std – standard deviation vector for d

Note

the system is modified to be solved by ordinary LS by the change G<- (G.T/std).T and d<- d/std

pyacs.lib.glinalg.lscov(G, d, cov, method='chol')[source]¶

Solve the least-squares (LS) problem with data covariance

Parameters

G – m x n model matrix as 2D numpy array
d – m 1D numpy observation vector
cov – covariance matrix for d

pyacs.lib.glinalg.lsw_full(G, d, std, verbose=False)[source]¶

Solve the least-squares (LS) with data uncertainties provided as a vector and returns the posterior covariance

Parameters

G – m x n model matrix as 2D numpy array
d – m 1D numpy observation vector
std – standard deviation vector for d

pyacs.lib.glinalg.lscov_full(G, d, cov, verbose=False)[source]¶

Solve the least-squares (LS) with data covariance and returns the posterior covariance

Parameters

G – m x n model matrix as 2D numpy array
d – m 1D numpy observation vector
cov – covariance matrix for d

pyacs.lib.glinalg.cov_to_invcov(M)[source]¶

Inverse a covariance matrix

Parameters: cov – 2D numpy array covariance matrix
Returns: 2D numpy array inverse covariance matrix

pyacs.lib.glinalg.corr_to_cov(corr, sigma_m)[source]¶

Correlation to covariance matrix

Parameters

corr – correlation matrix
sigma_m – vector of standard deviation = sqrt(diag(Cov))

:return covariance matrix

pyacs.lib.glinalg.cov_to_corr(Cov)[source]¶

Covariance to correlation transformation

Parameters: cov – covariance matrix
Return corr,sigma_m: correlation matrix and standard deviation matrix

pyacs.lib.glinalg.symmetrize(a, type_matrix)[source]¶

Extract the upper or lower triangle and make a symmetric matrix

Parameters

a – numpy 2D array, must be a square matrix
type – ‘triu’ or ‘tril’, triangle used to form the matrix

pyacs.lib.glinalg.dot_and_sum(LX, a)[source]¶: does a matrix by scalar product and sum it :param LX : list of arrays with the same dimension :param a : list of scalars

pyacs.lib.glinalg.repeat_matrix_in_col(G, n)[source]¶

Repeat a matrix along a column

R = np.empty( ( n,G.shape[0], G.shape[1] ) ) R[:] = G return( R.reshape( n*G.shape[0], G.shape[1]) )

pyacs.lib.glinalg.odot(a, G)[source]¶

Customize vector by matrix product.

Let’s assume we have a matrix made of equal dimension sub-matrices $G_1,cdots, G_n$ egin{equation} left[ egin{array}{c} G_1G_2

dots

G_n end{array}

ight]

end{equation} and a vector of scalars egin{equation} left[ egin{array}{c} a_1a_2

dots

a_n end{array}

ight]

end{equation} and we want egin{equation} left[ egin{array}{c} a_1 G_1a_2 G_2

dots

a_n G_n end{array}

ight]

end{equation} We do this with numpy broadcasting

param a: 1D numpy array of multipliers (scalars)
param G: 2-D numpy arrays of submatrix G_i.
return: result of the multiplication
note: if G.shape = (n,m), and a.shape=(l), then the shape of the G_i is (n/l,m)

pyacs.lib.glinalg.dot(A, B)[source]¶: Matrix/matrix, matrix/vector and vector/vector multiplication :author : P. Rebischung :date:Created : 02-Aug-2013 :changes : :param A: Matrix or vector :param B : Matrix or vector :return: A*B

pyacs.lib.glinalg.syminv(M)[source]¶: Invert a symmetric matrix :author : P. Rebischung :created : 02-Aug-2013 :param M : Matrix :return:Inverse matrix

pyacs.lib.glinalg.sympinv(M, verbose=False)[source]¶: Pseudo-invert a positive semi-definite symmetric matrix :author: P. Rebischung :created : 02-Aug-2013 :param M : Matrix :return: Pseudo-inverse

pyacs.lib.glinalg.make_normal_system(A, d, inv_Cd)[source]¶: Given the linear system A x = d with Cd the covariance matrix of d the associated normal system is A.T Cd-1 A x = A.T Cd-1 d. returns N= A.T Cd-1 A, Nd=A.T Cd-1 d

pyacs.lib.glinalg.make_normal_system_tarantola(G, d, m0, inv_Cd, inv_Cm)[source]¶: returns Gt inv_Cd G + inv_Cm , Gt inv_Cd d + inv_Cm m0

pyacs.lib.glinalg.matrix_from_pattern(pattern, structure)[source]¶

creates a matrix made of a pattern according to a structure

:param pattern : pattern 2D numpy array to be duplicated :param structure: 2D numpy array giving the structure as multiplicating factors :return: the block matrix as a 2D numpy array

>>> pattern = np.arange(6).reshape(2,3)
>>> structure = np.arange(6).reshape(3,2)

>>> print(pattern)
>>> print('--')
>>> print(structure)
>>> print('--')
>>> print( matrix_from_pattern( pattern , structure ) )

[[0 1 2] [3 4 5]] – [[0 1]

[2 3] [4 5]]

– [[ 0 0 0 0 1 2]

[ 0 0 0 3 4 5] [ 0 2 4 0 3 6] [ 6 8 10 9 12 15] [ 0 4 8 0 5 10] [12 16 20 15 20 25]]

pyacs.lib.glinalg.extract_block_diag(a, n, k=0)[source]¶

pyacs.lib.gmtpoint module¶

class pyacs.lib.gmtpoint.GMT_Point(code=None, lon=None, lat=None, he=0, Ve=None, Vn=None, Vu=0, SVe=0, SVn=0, SVu=0, SVen=0, Cv_xyz=None, Cv_enu=None, index=None)[source]¶

Bases: object

Simple Point Class

Attributes (mandatory) lon: longitude (decimal degrees) lat: latitude (decimal degrees)

Attributes (optinal): code : 4-letters code he : height above the ellipsoid ve, vn: east & north components of velocity sve,svn,sven: formal error (standard deviation) on velocity component and correlation [-1,1]

copy()[source]¶

get_info(display=False, legend=False, verbose=False)[source]¶: Print information of the current GMT_Point

magaz()[source]¶

Returns norm & azimuth of velocity from the current GMT_Point

Returns: norm,azimuth in mm/yr and decimal degree clockwise from North
Note: sigmas not implemented yet

assign_index(index)[source]¶

Assigns an index to the current GMT_Point. Useful for building linear systems.

Parameters: index – index

:return : the modified GMT_Point instance

get_index()[source]¶

Returns index of the current GMT_Point

Returns: index

spherical_distance(M)[source]¶

Returns spherical distance between two GMT_Point (in meters)

Parameters: M – GMT_Point instance

:return : distance along the sphere between current instance and M

pole(W=None, SW=None, type_euler='rot', option='predict')[source]¶: Substract the prediction of a given Euler pole

add_to_gmt_psvelo(fname, overwrite=False, verbose=False)[source]¶

Adds the current GMT_Point to a gmt psvelo file

Parameters

fname – gmt psvelo file
overwrite – will overwrite the line corresponding having the same code (default is False), generate an error otherwises
verbose – verbose mode

rotate_vel(angle, unit='radians')[source]¶

Rotate a vector of angle (clockwise)

Parameters

angle – angle
unit – ‘dec_deg’ or ‘radians’ (default radians)

:return : rotated components of velocity

midpoint(point, code='Mm')[source]¶

returns the mid-point between two GMT_Point

Parameters: point – GMT_Point

:return : mid_point as GMT_Point

pyacs.lib.gpstime module¶

Note: Leap seconds still needs to be checked an improve A Python implementation of GPS related time conversions.

Modifications for GPS seconds by Duncan Brown

PyUTCFromGpsSeconds added by Ben Johnson

This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

GPS Time Utility functions

This file contains a Python implementation of GPS related time conversions.

The two main functions convert between UTC and GPS time (GPS-week, time of week in seconds, GPS-day, time of day in seconds). The other functions are convenience wrappers around these base functions.

A good reference for GPS time issues is: http://www.oc.nps.navy.mil/~jclynch/timsys.html

Note that python time types are represented in seconds since (a platform dependent Python) Epoch. This makes implementation quite straight forward as compared to some algorigthms found in the literature and on the web.

pyacs.lib.gpstime.dayOfWeek(year, month, day)[source]¶: returns day of week: 0=Sun, 1=Mon, .., 6=Sat

pyacs.lib.gpstime.gpsWeek(year, month, day)[source]¶: returns (full) gpsWeek for given date (in UTC)

pyacs.lib.gpstime.julianDay(year, month, day)[source]¶: returns julian day=day since Jan 1 of year

pyacs.lib.gpstime.mkUTC(year, month, day, hour, minute, sec)[source]¶: similar to python’s mktime but for utc

pyacs.lib.gpstime.ymdhmsFromPyUTC(pyUTC)[source]¶: returns tuple from a python time value in UTC

pyacs.lib.gpstime.wtFromUTCpy(pyUTC, leapSecs=14)[source]¶: convenience function: allows to use python UTC times and returns only week and tow

pyacs.lib.gpstime.gpsFromUTC(year, month, day, hour, minute, ssec, leapSecs=30)[source]¶

converts UTC to: gpsWeek, secsOfWeek, gpsDay, secsOfDay

a good reference is: http://www.oc.nps.navy.mil/~jclynch/timsys.html

This is based on the following facts (see reference above):

GPS time is basically measured in (atomic) seconds since January 6, 1980, 00:00:00.0 (the GPS Epoch)

The GPS week starts on Saturday midnight (Sunday morning), and runs for 604800 seconds.

Currently, GPS time is 13 seconds ahead of UTC (see above reference). While GPS SVs transmit this difference and the date when another leap second takes effect, the use of leap seconds cannot be predicted. This routine is precise until the next leap second is introduced and has to be updated after that.

SOW = Seconds of Week SOD = Seconds of Day

Note: Python represents time in integer seconds, fractions are lost!!!

pyacs.lib.gpstime.UTCFromGps(gpsWeek, SOW, leapSecs=14)[source]¶

converts gps week and seconds to UTC

see comments of inverse function!

SOW = seconds of week gpsWeek is the full number (not modulo 1024)

pyacs.lib.gpstime.GpsSecondsFromPyUTC(pyUTC, _leapSecs=14)[source]¶

converts the python epoch to gps seconds

pyEpoch = the python epoch from time.time()

pyacs.lib.gpstime.testTimeStuff()[source]¶

pyacs.lib.gpstime.testJulD()[source]¶

pyacs.lib.gpstime.testGpsWeek()[source]¶

pyacs.lib.gpstime.testDayOfWeek()[source]¶

pyacs.lib.gpstime.testPyUtilties()[source]¶

pyacs.lib.icosahedron module¶

pyacs.lib.icosahedron.distance(A, B)[source]¶

pyacs.lib.icosahedron.icosahedron()[source]¶: Constructs an icosahedron on the unit-sphere

pyacs.lib.icosahedron.subdivide(verts, faces)[source]¶: Subdivide each triangle into four triangles, pushing verts to the unit sphere

pyacs.lib.icosahedron.mesh_global(num_subdivisions=6)[source]¶

pyacs.lib.icosahedron.mesh_regional(num_subdivisions=6, bounds=None)[source]¶: Makes an equilateral triangles mesh over a sphere of unit radius using successive subdivisions of an initial icosahedron returns verts and faces verts: list of vertices ; a vertice is a list of [x,y,z] in geocentric cartesian coordinates over the sphere of unit radius faces: list of faces ; a face is a list of vertices index faces[i][j] gives the [X, Y, Z] coordinates of vertex j of face i

pyacs.lib.pygamit_module module¶

module for Gamit

pyacs.lib.pygamit_module.read_qfile(qfile)[source]¶: Reads a qfile and returns the following information: - number of sites included in the analysis - list of sites that do not have solution - list of sites with large adjustment - total number of ambiguities - total number of fixed ambiguities - the pre/postfit values - Normal end in solve

pyacs.lib.pygamit_module.get_nearest_from_apr(aprfile, n, linclude, lexclude, X, Y, Z, dec_year)[source]¶: Returns the n nearest site from X,Y,Z present in linclude and apr

pyacs.lib.pygamit_module.substitute_site(aprfile, site, X, Y, Z, dec_year)[source]¶: Substitute an entry in apr file

pyacs.lib.pygamit_module.update_apr(input_apr, target_apr, threshold, action)[source]¶: Update an apr file using another apr threshold is a minimum value for updating if action is False only check

pyacs.lib.pygamit_module.make_apr_doy(input_apr, output_apr, lsite, doy, year)[source]¶: Creates an apr file for a given list of sites and select the best coordinates according to dates

pyacs.lib.robustestimators module¶

RobustEstimators.py includes a number of robust estimators to solve linear problems.

exception pyacs.lib.robustestimators.Error[source]¶

Bases: Exception

Base class for exceptions in module robustestimators.py

exception pyacs.lib.robustestimators.UnboundedFunctionError[source]¶

Bases: Exception

Exception raised for unbounded objective function.

pyacs.lib.robustestimators.Dikin(A, y, W, eps=0.003)[source]¶

L1-norm estimation of parameters x using the Dikin’s method using Linear optimization in linear model y=Ax+e

Parameters

A – Model matrix
y – observed values (observation vector)
W – Weight matrix of observables (of DIAGONAL type)
eps – small value for iteration (default: eps = 1e-6)

:return : x,e: vector of estimated parameters and residuals

Note: translated from Matlab code kindly provided by Amir Khodabandeh june.2009

reference:Recursive Algorithm for L1 Norm Estimation in Linear Models, A. Khodabandeh and A. R. Amiri-Simkooei, JOURNAL OF SURVEYING ENGINEERING ASCE / FEBRUARY 2011 / 1 doi:10.1061/ASCESU.1943-5428.0000031 Translated to python from Matlab original by J.-M. Nocquet July 2011

pyacs.lib.robustestimators.Dik_m(c, A, b, er)[source]¶: subject:solve the standard form of linear programming by affine/Dikin’s method(“an interior point method”) minimize z=c’*x; subject to Ax=b; input:(c):coefficients of objective function(z) as n-vector(hint:a column vector) (A):matrix of constraint set with size m*n (b):m-vector of constraint set (er):maximum discrepancy between two iteration.(“stopping criterion”) output:(X):unknown variables (z):optimal value of objective function (D):Centering transformer “D”(a diagonal matrix) Amir khodabandeh Oct.2008

pyacs.lib.shapefile module¶

Various subroutine to convert pyacs results to shapefiles

pyacs.lib.shapefile.static_slip_to_shapefile(static_slip_file, shp_name, dis_type='rec', verbose=False)[source]¶

Converts pyacs/pyeq slip solution into a shapefile

Parameters

static_slip_file – output from pyeqstatic_inversion.py (*sol_slip.dat)
shp_name – output shapefile name
dis_type – either ‘rec’ (rectangular dislocation) or ‘tde’ (triangular dislocation element)
verbose – verbose mode

pyacs.lib.shapefile.psvelo_to_shapefile(psvelo_file, shp_name, verbose=False)[source]¶

Converts a psvelo GMT file into shapefile

Parameters

psvelo – GMT psvelofile
shp_name – output shapefile name
verbose – verbose mode

pyacs.lib.shapefile.pyeblock_fault(pyeblock_file, shp_name, verbose=False)[source]¶

Converts a pyeblock fault file into a shapefile

Parameters

pyeblock_fault – pyeblock fault file
shp_name – output shapefile name
verbose – verbose mode

pyacs.lib.strain module¶

Strain rate calculation library

pyacs.lib.strain.vgrad(X, Y, VE, VN, SVE=None, SVN=None, CVEN=None, CODE=None, method='WLS', verbose=False)[source]¶

Linear estimates of a horizontal velocity gradient

Parameters

X,Y – 1D numpy array of longitudes and latitudes
VE,VN – 1D numpy array of east and north velocity. Unit is assummed to be mm/yr
SVE,SVN,CVEN – 1D numpy array of east and north velocity standard deviation (mm/yr) and correlation.
CODE – 1D numpy of string with site codes
method – estimator in ‘L1’,’WLS’ (default: weighted least-squares ‘WLS’)

Returns

DV, the velocity gradient tensor as a 2x2 2D-numpy array and a 2-columns (East and North) residuals 2D numpy array

pyacs.lib.syscmd module¶

Runs a system commands

exception pyacs.lib.syscmd.Error[source]¶

Bases: Exception

Base class for exceptions in module syscmd.

exception pyacs.lib.syscmd.ReturnedCodeError[source]¶

Bases: Exception

Exception raised for returned code not 0.

pyacs.lib.syscmd.getstatusoutput(cmd, verbose=False)[source]¶

pyacs.lib.syscmd.run_command_verbose(command)[source]¶

pyacs.lib.timeperiod module¶

exception pyacs.lib.timeperiod.TimePeriodUndefinedError[source]¶: Bases: Exception

class pyacs.lib.timeperiod.TimePeriod(start_time=None, end_time=None)[source]¶

Bases: object

The beginning of this period as datetime object

isdefined()[source]¶

begin()[source]¶

epoch_begin()[source]¶

end()[source]¶

epoch_end()[source]¶

intersection(period)[source]¶

has_in(date)[source]¶

Tells whether a given date is in the period

return: True if the given date in within the period, False otherwise

display()[source]¶

get_info()[source]¶

pyacs.lib.units module¶

Various routines dealing with units conversion.

pyacs.lib.units.mas2rad(mas)[source]¶

Converts milliarcseconds to radians

Parameters: mas – angle in millarcseconds
Returns: rad: in radians
Example

>>> import pyacs
>>> pyacs.mas2rad(1.)
>>> 4.84813681109536e-09
>>> pyacs.mas2rad(np.array([1.0,2.0]))
>>> array([  4.84813681e-09,   9.69627362e-09])

pyacs.lib.units.rad2mas(rad)[source]¶

Converts radians to milliarcseconds

Parameters: rad – angle in radians
Returns mas: angle in millarcseconds:

pyacs.lib.units.radians2deg_mn_sec(rad, angle_range='-180-180')[source]¶

Converts an angle from radians to deg, minutes, seconds

Parameters

rad – angle in radians
angle_range – either ‘-180-180’ or ‘0-360’

Returns

deg, minutes, seconds

Return type

int,int,float

pyacs.lib.units.moment_to_magnitude(moment)[source]¶

Converts a moment in N.m to Mw

Parameters: moment – in N.m
Return magnitude: 2./3.*(math.log10(M0)-9.05)

pyacs.lib.units.magnitude_to_moment(magnitude)[source]¶

Converts a Mw magnitude to moment in N.m

Parameters: magnitude – Mw
Return moment: in N.m, 10^(1.5*magnitude+9.5)

pyacs.lib.utils module¶

Various useful routines

pyacs.lib.utils.numpy_array_2_numpy_recarray(A, names)[source]¶: Converts a numpy array to a numpy recarray names is the names of each field

pyacs.lib.utils.numpy_recarray_2_numpy_array(A)[source]¶: Converts a structured array (with homogeneous dtype) to a np.array

pyacs.lib.utils.save_np_array_with_string(A, S, fmt, outfile, comment='')[source]¶

saves a numpy array as a text file. This is equivalent to np.savetxt except that S contains strings that are added to each row.

Parameters

A – 2D numpy array to be saved
S – 1D numpy string array. S.ndim = A.shape[0]
fmt – format for printing
outfile – out file

Example

A = np.arange(4).reshape(-1,2) S = np.array([‘lima’,’quito’]) from pyacs.lib.utils import save_np_array_with_string save_np_array_with_string(A,S,”%03d %03d %s”,’test.dat’)

pyacs.lib.utils.make_grid(min_lon, max_lon, min_lat, max_lat, nx=None, ny=None, step_x=None, step_y=None, outfile=None, format='psvelo', comment='')[source]¶

Generates a text file as a grid

Parameters

min_lon,max_lon,min_lat,max_lat – grid bounds coordinates in decimal degrees
ny (nx) – number of points in the grid along longitude. If ny is not provided, ny=nx
step_x,step_y – step for the grid. This is an alternative to nx. If step_y is not provided, step_y=step_x
outfile – output file name. Default=None, no file written
format – if format is None, then only lon, lat are written. If format=’psvelo’ (default), then the line is filled with 0. and sequentially site names
comment – comment to be added to the output file.

Returns

the grid as 2D numpy array

pyacs.lib.utils.str2list_float(my_str)[source]¶

Converts a list provided as a string to a true list of float

Parameters: my_str – string e.g. ‘[0,2.5,1E9]’

:return:[0,2.5,1E9]

pyacs.lib.vectors module¶

Useful operations on vectors (1-D numpy array)

pyacs.lib.vectors.np_array(X)[source]¶

Converts to 1-D numpy array

Parameters: X – list or tuple to be converted

:return : 1D numpy array

pyacs.lib.vectors.norm(X)[source]¶

Returns the norm of a vector

Parameters: X – list, tuple or 1D numpy array

:return : norm (float)

pyacs.lib.vectors.normalize(X)[source]¶

Normalize a vector

Parameters: X – list, tuple or 1D numpy array

:return : 1D numpy array

pyacs.lib.vectors.scal_prod(X, Y)[source]¶

Returns the scalar products

Parameters: X,Y – list, tuple or 1D numpy array

:return : scalar product as float

pyacs.lib.vectors.vector_product(A, B)[source]¶

Vector product

A[1]*B[2] - B[1]*A[2]

vector_product(A,B)= A[2]*B[0] - B[2]*A[0]: A[0]*B[1] - B[0]*A[1]

Parameters: A,B – 1D numpy arrays
Return AxB: 1D numpy array

pyacs.lib.vel_field module¶

class pyacs.lib.vel_field.Velocity_Field(file_name=None, name=None, lgmt_points=None, verbose=False)[source]¶

Bases: object

Velocity_Field class: reads a velocity field from a GMT psvelo file and provides methods to manipulate velocity field

classmethod read(file_name=None, lexclude=[], lonly=[], verbose=False)[source]¶: Reads a GMT psvelo file

info(details=True)[source]¶: Print basics information

add_point(M)[source]¶: Appends a GMT_Point to a velocity field

remove_point(code)[source]¶

Removes the a GMT_Point to a velocity field given its code

Parameters: code – 4-characters code for the point to be removed

write(out_file=None, lexclude=[], verbose=True, comment='', up=False)[source]¶: Writes a GMT psvelo file a list site name to be excluded can be provided

nsites()[source]¶: Returns the number of sites in velocity field

l_GMT_Point()[source]¶: Returns the velocity field as a list of GMT_Point object

print_info_site(code, verbose=False)[source]¶: Print the info for a given site by his code

subset(lonly=None, lexclude=None)[source]¶

Returns a new Velocity_Field object from a subset of sites

Parameters: lonly,lexclude – list of site codes to be included or excluded

radial(center)[source]¶

Returns a Velocity Field whose components are radial and tangential with respect to a given center

Parameters: center – numpy array [longitude,latitude] of center

lcode()[source]¶: Returns a list of all point codes in current Velocity Field object

site(code)[source]¶: Returns a site in current Field object as a GMT_Point object

calc_pole(plates, method='WLS', verbose=False)[source]¶

Performs Euler pole calculation

Parameters

plates – dictionary with the name of plate as key and list of sites for each plate as values
method – choose among weighted least-squares ‘WLS’ or ‘L1’

Output

for each plate, the following files will be created: - euler_stat_plate_name.dat: Euler pole values and statistics - plate_name.gmt: velocities (GMT psvelo format) with respect to the calculated Euler pole - plate_name.shp: velocities (shapefile format) with respect to the calculated Euler pole - eulers_sum.dat: summary of Euler pole values

pole(lexclude=[], method='WLS', exp='pLate', log=None)[source]¶: Euler pole calculation; available optimization methods are WLS: weighted least squares, LS: least-squares, Dikin: L1 linear estimator

substract_pole(W=None, type_euler='rot')[source]¶: substract the prediction of an Euler pole to the velocity field

common(linGpoint, prefit=10.0, lexclude=[])[source]¶: Returns a list of sites common to the current Velocity Field object and the list of GMT_Points provided in argument Coordinates will be the ones from the Sinex and NOT from the list of Gpoints Commons point are points with same code pt and soln

proj_profile(slon, slat, elon, elat, d, save=None, verbose=False)[source]¶

project velocity components along a great circle defined by initial/stop points (slon,slat,elon,elat)

Parameters

slon,slat – coordinates of profile start point (decimal degrees)
elon,elat – coordinates of profile end point (decimal degrees)

:param d : maximum distance for a point to be considered :param save : output file name (optional)

:return : numpy 1D array with np_code, np_distance_along_profile, np_distance_to_profile , np_Ve , np_Vn , np_SVe , np_SVn , np_v_parallele , np_v_perpendicular , np_sigma_v_parallele , np_sigma_v_perpendicular , np_lazimuth

strain(lcode, save=None, method='WLS', verbose=False)[source]¶

Calculates the strain rate given a list of sites

Parameters: lcode – list of site names

:param save : file to save the result :param method: estimator in ‘L1’,’WLS’ (default: weighted least-squares ‘WLS’) :param verbose: verbose mode

pyacs.lib package¶

Subpackages¶

Submodules¶

pyacs.lib.astrotime module¶

pyacs.lib.coordinates module¶

pyacs.lib.errors module¶

pyacs.lib.euclid module¶

pyacs.lib.euler module¶

pyacs.lib.glinalg module¶

pyacs.lib.gmtpoint module¶

pyacs.lib.gpstime module¶

pyacs.lib.icosahedron module¶

pyacs.lib.pygamit_module module¶

pyacs.lib.robustestimators module¶

pyacs.lib.shapefile module¶

pyacs.lib.strain module¶

pyacs.lib.syscmd module¶

pyacs.lib.timeperiod module¶

pyacs.lib.units module¶

pyacs.lib.utils module¶

pyacs.lib.vectors module¶

pyacs.lib.vel_field module¶

Module contents¶