[GRASSLIST:4183] Re: structural geology - OFF TOPIC
Phillip J. Allen
paallen at attglobal.net
Sun Aug 15 09:24:02 EDT 2004
Carlos,
I have written a function in pl/pgSQL (PostgreSQL) and Visual Basic (for
Access) that takes Dip-Azimuth (0- <360 horizontal) and Vertical Dip
Angle (0-90 downwards) and recalculates coordinates for a stereonet.
The sterionet coordinates have the origin 0,0 and a circle radius of
100 unites. So with your structural/orientation data in a database use
this function in a query and open the query in any graphing program and
you have a stereonet. I am not really sure if it is a wolff or schmitt
net but I could never really tell them appart anyway. I always kept
meaning to work these functions up into a R function.
I have drawn the sterionet lines in Autocad and just open the dwg up in
ArcView and can then sybolize the poles with any attribute I need.
Below are the PostgreSQL and Access functions. Good luck.
I have always been meaning to work up a little script in
PostgreSQL/PostGIS to create a 3D-Circle or Thin Disk representing the
structural point and it's orientation ot view in 3-d but I have just
never gotten around to it yet.
Good Luck.
Phillip J. Allen
Chief Geochemist Anglogoldashanti South America
Lima Peru
e-mail: paallen at attglobal.net
Carlos Henrique Grohmann wrote:
>Hello, sorry about the off-topic, but dows anyone kkow a software for structural
>geology (stereonet) for linux? I have beeen looking on the net, but haven't
>find nothing yet.
>
>thanks
>
>Carlos
>
>
****************************
*****************************
**PostgreSQL PL/PGSQL Function**
*****************************
*****************************
-- Name: "stereonet" (smallint,double precision,double precision) Type:
FUNCTION Owner: paallen
--
CREATE FUNCTION "stereonet" (smallint,double precision,double precision)
RETURNS double precision AS '
DECLARE
xyz ALIAS FOR $1;
dipazimuth FLOAT8;
dipangle FLOAT8;
originx CONSTANT FLOAT8 := 0.0;
originy CONSTANT FLOAT8 := 0.0;
circdiam CONSTANT FLOAT8 := 100.0;
spherex FLOAT8;
spherey FLOAT8;
spherez FLOAT8;
stnetx FLOAT8;
stnety FLOAT8;
BEGIN
-- Calculate XY Coordinates for Wulff Steronet with origin at 0,0,0 and
circle diameter of 100
-- Convert Plane to Pole
dipazimuth = $2;
dipangle = $3;
dipazimuth = dipazimuth + 180.0;
IF dipazimuth >= 360.0 THEN
dipazimuth = dipazimuth -360.0;
END IF;
dipangle = 0 - (90.0 - dipangle);
spherex = returnpc2newxyz(originx,
originy,0.0,circdiam,dipazimuth, dipangle,1);
spherey = returnpc2newxyz(originx,
originy,0.0,circdiam,dipazimuth,dipangle,2);
spherez = returnpc2newxyz(originx,
originy,0.0,circdiam,dipazimuth,dipangle,3);
IF spherez >=0 THEN
stnetx = spherex;
stnety = spherey;
ELSE
stnetx = (((0 - circdiam) / (spherez - 100)) * (spherex -
originx)) + originx;
stnety = (((0 - circdiam) / (spherez - 100)) * (spherey -
originy)) + originy;
END IF;
IF xyz = 1 THEN
RETURN stnetx;
ELSIF xyz = 2 THEN
RETURN stnety;
ELSIF xyz = 3 THEN
RETURN 0.0;
ELSIF xyz = 4 THEN
RETURN spherex;
ELSIF xyz = 5 THEN
RETURN spherey;
ELSIF xyz = 6 THEN
RETURN spherez;
ELSE
RETURN -9999999.0;
END IF;
END;
' LANGUAGE 'plpgsql';
--
-- Name: "returnpc2newxyz" (double precision,double precision,double
precision,double precision,double precision,double precision,smallint)
Type: FUNCTION Owner: paallen
--
CREATE FUNCTION "returnpc2newxyz" (double precision,double
precision,double precision,double precision,double precision,double
precision,smallint) RETURNS double precision AS '
DECLARE
newx FLOAT8;
newy FLOAT8;
newz FLOAT8;
xcenter ALIAS FOR $1;
ycenter ALIAS FOR $2;
zcenter ALIAS FOR $3;
dist ALIAS FOR $4;
azimuthdd ALIAS FOR $5;
inclinangdd ALIAS FOR $6;
horizdist FLOAT8;
ddrad CONSTANT FLOAT8 := 3.14159265358979/180.0;
BEGIN
IF inclinangdd < 0 THEN
newz = zcenter - (dist * sin((abs(inclinangdd)) * ddrad));
ELSE
newz = zcenter + (dist * sin((abs(inclinangdd)) * ddrad));
END IF;
IF $7 = 3 THEN
RETURN newz;
END IF;
horizdist = dist * cos(abs(inclinangdd * ddrad));
IF $7 = 1 THEN
newx = xcenter + (horizdist * (sin(azimuthdd * ddrad)));
RETURN newx;
ELSIF $7 = 2 THEN
newy = ycenter + (horizdist * (cos(azimuthdd * ddrad)));
RETURN newy;
END IF;
RETURN -9999999.0;
END;
' LANGUAGE 'plpgsql';
*****************************
*****************************
**Access Visual Basic Functions**
*****************************
*****************************
Function StereoNet(ByVal snt As String, xyz As String, dazimuth As
Double, dangle As Double) As Double
Rem*************************************************************************************************************
Rem*********Calculate Wulff or Lambert X Stereo Net Coordinates from Dip
Azimuth & Dip Angle degrees************
Rem*************************************************************************************************************
Rem****************************Great Circle is 100 units
Radius*************************************************
Dim originx, originy, circdiam As Double
Dim spherex, spherey, spherez As Double
Dim stnetx, stnety As Double
Dim anb As Integer
originx = 0
originy = 0
circdiam = 100
Rem Convert plane to pole
dazimuth = dazimuth + 180
If (dazimuth >= 360) Then
dazimuth = dazimuth - 360
End If
dangle = 90 - dangle
Rem anb = MsgBox(dazimuth, vbOKOnly, "azimuth")
Rem anb = MsgBox(dangle, vbOKOnly, "angle")
Rem If (snt = "W") Then Rem Wulff Stereo Net
spherex = ReturnPC2NewX(originx, originy, 0, circdiam, dazimuth, 0 -
dangle)
spherey = ReturnPC2NewY(originx, originy, 0, circdiam, dazimuth, 0 -
dangle)
spherez = ReturnPC2NewZ(originx, originy, 0, circdiam, dazimuth, 0 -
dangle)
Rem anb = MsgBox(spherex, vbOKOnly, "spherex")
Rem anb = MsgBox(spherey, vbOKOnly, "spherey")
Rem anb = MsgBox(spherez, vbOKOnly, "spherez")
If (spherez >= 0) Then
stnetx = spherex
stnety = spherey
Else
stnetx = (((0 - circdiam) / (spherez - 100)) * (spherex -
originx)) + originx
stnety = (((0 - circdiam) / (spherez - 100)) * (spherey -
originy)) + originy
End If
Rem ElseIf (snt = "L") Then Rem Lambert Stereo Netdssd
Rem End If
If xyz = "x" Then
StereoNet = stnetx
ElseIf xyz = "y" Then
StereoNet = stnety
ElseIf xyz = "sx" Then
StereoNet = spherex
ElseIf xyz = "sy" Then
StereoNet = spherey
ElseIf xyz = "sz" Then
StereoNet = spherez
End If
End Function
Public Function ReturnPC2NewY(ByVal XCPoint, YCPoint, ZCPoint, dist,
AzimuthAnglDD, InclinAnglDD)
Dim NewZ As Double
Dim DH As Double
Rem Calculate Depth(z)+z is upward, 0 is horizontal, -z is downward
If InclinAnglDD < 0 Then
NewZ = ZCPoint - (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979
/ 180))
Else
NewZ = ZCPoint + (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979
/ 180))
End If
Rem Calculate Horizontal X Distance
DH = dist * Cos(Abs(InclinAnglDD * 3.14159265358979 / 180))
ReturnPC2NewY = YCPoint + (DH * (Cos(AzimuthAnglDD * 3.14159265358979 /
180)))
End Function
Function ReturnPC2NewZ(ByVal XCPoint, YCPoint, ZCPoint, dist,
AzimuthAnglDD, InclinAnglDD)
Dim DH As Double
Rem Calculate Depth(z)+z is upward, 0 is horizontal, -z is downward
If InclinAnglDD < 0 Then
ReturnPC2NewZ = ZCPoint - (dist * Sin((Abs(InclinAnglDD)) *
3.14159265358979 / 180))
Else
ReturnPC2NewZ = ZCPoint + (dist * Sin((Abs(InclinAnglDD)) *
3.14159265358979 / 180))
End If
End Function
Function ReturnPC2NewX(ByVal XCPoint, YCPoint, ZCPoint, dist,
AzimuthAnglDD, InclinAnglDD)
Dim NewX As Double, NewZ As Double
Dim DH As Double
Rem Calculate Depth(z)+z is upward, 0 is horizontal, -z is downward
If InclinAnglDD < 0 Then
NewZ = ZCPoint - (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979
/ 180))
Else
NewZ = ZCPoint + (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979
/ 180))
End If
Rem Calculate Horizontal X Distance
DH = dist * Cos(Abs(InclinAnglDD * 3.14159265358979 / 180))
ReturnPC2NewX = XCPoint + (DH * (Sin(AzimuthAnglDD * 3.14159265358979 /
180)))
End Function
More information about the grass-user
mailing list