/[escript]/trunk/modellib/py_src/crustal/mines.py
ViewVC logotype

Annotation of /trunk/modellib/py_src/crustal/mines.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 981 - (hide annotations)
Mon Feb 19 06:52:25 2007 UTC (14 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 19903 byte(s)
extraction records added
1 gross 977 """
2     mining activities
3    
4     @var __author__: name of author
5     @var __licence__: licence agreement
6     @var __url__: url entry point on documentation
7     @var __version__: version
8     @var __date__: date of the version
9     """
10    
11     __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
12     http://www.access.edu.au
13     Primary Business: Queensland, Australia"""
14     __license__="""Licensed under the Open Software License version 3.0
15     http://www.opensource.org/licenses/osl-3.0.php"""
16     __author__="Lutz Gross, l.gross@uq.edu.au"
17     __url__="http://www.iservo.edu.au/esys/escript"
18     __version__="$Revision$"
19     __date__="$Date$"
20    
21     from xml.dom import minidom
22     from esys.pycad import *
23     from math import pi, cos, sin, sqrt
24     from esys.pycad.gmsh import Design
25     #======== that should be in a different file =============================
26     class LocationOnEarth(object):
27     AGD84="AGD84"
28     AGD84_semi_major_axis=6378160.00
29     AGD84_inverse_flatening=298.2500
30     DEG_UNIT="D.MM"
31     LENGTH_UNIT="M"
32     def __init__(self, longitude=0., latitude=0., altitude=0.):
33     """
34     """
35     self.__longitude=longitude
36     self.__latitude=latitude
37     self.__altitude=altitude
38     def getLLA(self):
39     return (self.__longitude,self.__latitude,self.__altitude)
40    
41     def getRelativeXYZ(self,center):
42     """
43     return the location in global coordinates relative to the center of the earth/
44     default reference grid is AGD84.
45    
46     @note: formula for AGD84 from http://www.ga.gov.au/geodesy/datums/transxyz.xls
47     """
48     C_LONG,C_LAT,C_ALT=center.getLLA()
49     C_LONG=unitConverter(C_LONG, center.DEG_UNIT,"RAD")
50     C_LAT=unitConverter(C_LAT, center.DEG_UNIT,"RAD")
51     I7=unitConverter(self.__latitude,self.DEG_UNIT,"RAD")
52     I8=unitConverter(self.__longitude,self.DEG_UNIT,"RAD")
53     C_ALT=unitConverter(C_ALT, center.LENGTH_UNIT,self.LENGTH_UNIT)
54     X=self.AGD84_semi_major_axis*sin(I7-C_LAT)
55     Y=self.AGD84_semi_major_axis*sin(I8-C_LONG)
56     Z=(self.__altitude-C_ALT)
57     return numarray.array([X,Y,Z])
58    
59    
60     def getXYZ(self,center=None,reference_grid=None):
61     """
62     return the location in global coordinates relative to the center of the earth/
63     default reference grid is AGD84.
64    
65     @note: formula for AGD84 from http://www.ga.gov.au/geodesy/datums/transxyz.xls
66     """
67     if reference_grid==None:
68     I7=unitConverter(self.__latitude,self.DEG_UNIT,"RAD")
69     I8=unitConverter(self.__longitude,self.DEG_UNIT,"RAD")
70     D15=1./self.AGD84_inverse_flatening
71     D16=2*D15-D15*D15
72     D17=self.AGD84_semi_major_axis/sqrt(1-D16*sin(I7)*sin(I7))
73     X=(D17 +self.__altitude)*cos(I7)*cos(I8)
74     Y=(D17 +self.__altitude)*cos(I7)*sin(I8)
75     Z=((1-D16)*D17+self.__altitude)*sin(I7)
76     else:
77     raise ValueError("reference grid %s is not supported."%reference_grid)
78     return numarray.array([X,Y,Z])
79    
80     def unitConverter(val,val_unit,out_unit):
81     VAL_UNIT=val_unit.upper()
82     OUT_UNIT=out_unit.upper()
83     out=None
84     if OUT_UNIT == "D.MM":
85     if VAL_UNIT == "D.MM":
86     out=float(val)
87     elif OUT_UNIT == "RAD":
88     if VAL_UNIT == "D.MM":
89     out=float(val)/180.*pi
90     elif OUT_UNIT == "M":
91     if VAL_UNIT == "M":
92     out=float(val)
93     elif VAL_UNIT == "KM":
94     out=float(val)*1000.
95 gross 981 elif OUT_UNIT == "KG":
96     if VAL_UNIT == "MT":
97     out=float(val)*1e9
98 gross 977 if out == None:
99     raise ValueError("cannot convert from unit %s to %s."%(val_unit,out_unit))
100     return out
101     #======== end of snip ====================================================
102     class BoxInTheCrust(object):
103     """
104     defines a box in the outer crust
105     """
106     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
107     self.__longitude_length=longitude_length
108     self.__latitude_length=latitude_length
109     self.__name=name
110     self.__location=location
111     self.__tilt=tilt
112     self.__depth=depth
113     self.__centerxyz=self.__location.getXYZ()
114     def getDiameter(self):
115     return sqrt(self.__depth**2+self.__longitude_length**2+self.__latitude_length**2)
116     def getCenterXYZ(self):
117     return self.__centerxyz
118     def getCenter(self):
119     return self.__location
120     def getName(self):
121     return self.__name
122     def getLongitudeLength(self):
123     return self.__longitude_length
124     def getLatitudeLength(self):
125     return self.__latitude_length
126     def getDepth(self):
127     return self.__depth
128     def getTotalDepth(self):
129     return self.__depth-self.__location.getLLA()[2]
130    
131     def createSurfaceLoop(self):
132     """
133     this creates a L{SurfaceLoop} describing the region
134     """
135     p1=Point(-self.__longitude_length/2,-self.__latitude_length/2,-self.__depth)
136     p2=Point( self.__longitude_length/2, self.__latitude_length/2,0.)
137     out=Brick(p1,p2)
138     out*=Rotatation(axis=[0,0,1], angle=self.__tilt)
139     return out
140    
141    
142     class MiningArea(BoxInTheCrust):
143     """
144     defines minuing activites in an area
145     """
146     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0., mines=[]):
147     super(MiningArea, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
148     self.__mines=mines
149 gross 981 def getStartOfRecords(self):
150     return min([ m.getStartOfRecords() for m in self.__mines ])
151 gross 977 def fillDesign(self,design):
152     """
153     this puts the mining area and all the mines into the given design
154     """
155     # define a bounding box around the mines:
156     X_MIN=0
157     X_MAX=0
158     Y_MIN=0
159     Y_MAX=0
160     DEPTH_MAX=0.
161     for m in self.__mines:
162     long=m.getLongitudeLength()
163     lat=m.getLatitudeLength()
164     pos=m.getCenter().getRelativeXYZ(self.getCenter())
165     X_MIN=min(X_MIN,pos[0]-lat/2.)
166     X_MAX=max(X_MAX,pos[0]+lat/2.)
167     Y_MIN=min(Y_MIN,pos[1]-long/2.)
168     Y_MAX=max(Y_MAX,pos[1]+long/2.)
169     DEPTH_MAX=max(DEPTH_MAX,m.getDepth()-pos[2])
170     lx=(X_MAX-X_MIN)/2*1.6
171     cx=(X_MIN+X_MAX)/2
172     X_MAX=cx+lx
173     X_MIN=cx-lx
174     ly=(Y_MAX-Y_MIN)/2*1.6
175     cy=(Y_MIN+Y_MAX)/2
176     Y_MAX=cy+ly
177     Y_MIN=cy-ly
178     DEPTH_MAX*=2.
179     dx=[X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX]
180     bb_p000=Point(X_MIN,Y_MIN,-DEPTH_MAX)
181     bb_p100=bb_p000+[dx[0],0.,0.]
182     bb_p010=bb_p000+[0.,dx[1],0.]
183     bb_p110=bb_p000+[dx[0],dx[1],0.]
184     bb_p001=bb_p000+[0.,0.,dx[2]]
185     bb_p101=bb_p000+[dx[0],0.,dx[2]]
186     bb_p011=bb_p000+[0.,dx[1],dx[2]]
187     bb_p111=bb_p000+[dx[0],dx[1],dx[2]]
188     bb_l10=Line(bb_p000,bb_p100)
189     bb_l20=Line(bb_p100,bb_p110)
190     bb_l30=Line(bb_p110,bb_p010)
191     bb_l40=Line(bb_p010,bb_p000)
192     bb_l11=Line(bb_p000,bb_p001)
193     bb_l21=Line(bb_p100,bb_p101)
194     bb_l31=Line(bb_p110,bb_p111)
195     bb_l41=Line(bb_p010,bb_p011)
196     bb_l12=Line(bb_p001,bb_p101)
197     bb_l22=Line(bb_p101,bb_p111)
198     bb_l32=Line(bb_p111,bb_p011)
199     bb_l42=Line(bb_p011,bb_p001)
200     bb_bottom=PlaneSurface(CurveLoop(-bb_l10,-bb_l40,-bb_l30,-bb_l20))
201     bb_top=PlaneSurface(CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42))
202     bb_front=PlaneSurface(CurveLoop(-bb_l11,bb_l10,bb_l21,-bb_l12))
203     bb_back=PlaneSurface(CurveLoop(bb_l30,bb_l41,-bb_l32,-bb_l31))
204     bb_left=PlaneSurface(CurveLoop(bb_l11,-bb_l42,-bb_l41,bb_l40))
205     bb_right=PlaneSurface(CurveLoop(-bb_l21,bb_l20,bb_l31,-bb_l22))
206     bb=SurfaceLoop(bb_bottom,bb_top,bb_front,bb_back,bb_left,bb_right)
207     bb.setLocalScale(min(X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX)/design.getElementSize())
208     # put all mines in to the domain:
209     mboxes=[]
210     props=[]
211     for m in self.__mines:
212     mb=m.createSurfaceLoop()
213     mb.setLocalScale(m.getDiameter()/design.getElementSize())
214     mb+=m.getCenter().getRelativeXYZ(self.getCenter())
215     mboxes.append(mb)
216     props.append(PropertySet(m.getName(),Volume(mb)))
217     # design.addItems(*tuple(props))
218     design.addItems(Volume(bb, holes= mboxes))
219     long=self.getLongitudeLength()
220     lat=self.getLatitudeLength()
221     dep=self.getDepth()
222     p000=Point(-long/2,-lat/2,-dep)
223     p100=p000+[long,0.,0.]
224     p010=p000+[0.,lat,0.]
225     p110=p000+[long,lat,0.]
226     p001=p000+[0.,0.,dep]
227     p101=p000+[long,0.,dep]
228     p011=p000+[0.,lat,dep]
229     p111=p000+[long,lat,dep]
230     l10=Line(p000,p100)
231     l20=Line(p100,p110)
232     l30=Line(p110,p010)
233     l40=Line(p010,p000)
234     l11=Line(p000,p001)
235     l21=Line(p100,p101)
236     l31=Line(p110,p111)
237     l41=Line(p010,p011)
238     l12=Line(p001,p101)
239     l22=Line(p101,p111)
240     l32=Line(p111,p011)
241     l42=Line(p011,p001)
242     bottom=PlaneSurface(CurveLoop(-l10,-l40,-l30,-l20))
243     top=PlaneSurface(CurveLoop(l12,l22,l32,l42),holes=[CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42)])
244     front=PlaneSurface(CurveLoop(-l11,l10,l21,-l12))
245     back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
246     left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
247     right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
248     dom=SurfaceLoop(bottom,top,front,back,left,right,-bb_bottom,-bb_front,-bb_back,-bb_left,-bb_right)
249     design.addItems(Volume(dom))
250     return
251    
252     class Mine(BoxInTheCrust):
253     """
254     defines a mine
255     """
256     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
257     super(Mine, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
258 gross 981 self.__record=[]
259 gross 977
260 gross 981 def addRecord(self, material, year, extraction):
261     TYPE=material.upper()
262     for y, e in self.__record:
263     if y == year:
264     if e.has_key(TYPE):
265     e[TYPE]+=extraction
266     else:
267     e[TYPE]=extraction
268     return
269     if y>year:
270     self.__record.insert(self.__record.index((y,e)), { TYPE : extraction} )
271     return
272     self.__record.append((year, { TYPE : extraction} ))
273     return
274     def getStartOfRecords(self):
275     if len(self.__record)>0:
276     return self.__record[0][0]
277     else:
278     raise ValueError("empty record of %s mine."%self.getName())
279 gross 977
280 gross 981
281 gross 977 def _parse(root):
282     """
283     parses child roots on root
284     """
285     if isinstance(root, minidom.Element):
286     if root.tagName == 'MiningArea':
287     NAME="unknown"
288     LOC=LocationOnEarth()
289     TILT=0.
290     LONGLENGTH=0.
291     LATLENGTH=0.
292     DEPTH=0.
293     MINES=[]
294     for node in root.childNodes:
295     if isinstance(node, minidom.Element):
296     if node.tagName == 'Tilt': TILT=_parse(node)
297     if node.tagName == 'Location': LOC=_parse(node)
298     if node.tagName == 'Name': NAME=_parse(node)
299     if node.tagName == 'Mine': MINES.append(_parse(node))
300     if node.tagName == 'Depth': DEPTH=_parse(node)
301     if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
302     if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
303     return MiningArea(location=LOC, name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, tilt=TILT, mines=MINES, depth=DEPTH)
304     elif root.tagName == 'Mine':
305     NAME="unknown"
306     LOC=LocationOnEarth()
307     TILT=0.
308     DEPTH=0.
309     LONGLENGTH=0.
310     LATLENGTH=0.
311 gross 981 RECORD=[]
312 gross 977 for node in root.childNodes:
313     if isinstance(node, minidom.Element):
314     if node.tagName == 'Name': NAME=_parse(node)
315     if node.tagName == 'Location': LOC=_parse(node)
316     if node.tagName == 'Tilt': TILT=_parse(node)
317     if node.tagName == 'Depth': DEPTH=_parse(node)
318     if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
319     if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
320 gross 981 if node.tagName == 'Record': RECORD.append(_parse(node))
321     m=Mine(location=LOC, tilt=TILT,name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, depth=DEPTH)
322     for r in RECORD:
323     m.addRecord(material=r[0],year=r[1],extraction=r[2])
324     return m
325 gross 977 elif root.tagName == 'LocationOnEarth':
326     long=0.
327     lat=0.
328     alt=0.
329     for node in root.childNodes:
330     if isinstance(node, minidom.Element):
331     if node.tagName == 'Longitude': long=_parse(node)
332     if node.tagName == 'Latitude': lat=_parse(node)
333     if node.tagName == 'Altitude': alt=_parse(node)
334     return LocationOnEarth(longitude=long, latitude=lat, altitude=alt)
335 gross 981 elif root.tagName == 'Record':
336     year=0
337     mat="unknown"
338     ext=0.
339     for node in root.childNodes:
340     if isinstance(node, minidom.Element):
341     if node.tagName == 'Year': year=_parse(node)
342     if node.tagName == 'Material': mat=_parse(node)
343     if node.tagName == 'Extraction': ext=_parse(node)
344     return (mat,year,ext)
345 gross 977 elif root.tagName == 'SouthWest':
346     return _parse(root.getElementsByTagName('LocationOnEarth')[0])
347     elif root.tagName == 'NorthEast':
348     return _parse(root.getElementsByTagName('LocationOnEarth')[0])
349     elif root.tagName == 'Longitude':
350     if root.hasAttribute("unit"):
351     unit=root.getAttribute("unit")
352     else:
353     unit="D.MM"
354     for node in root.childNodes:
355     if isinstance(node, minidom.Text):
356     return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
357     return 0.
358     elif root.tagName== 'Latitude':
359     if root.hasAttribute("unit"):
360     unit=root.getAttribute("unit")
361     else:
362     unit="D.MM"
363     for node in root.childNodes:
364     if isinstance(node, minidom.Text):
365     return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
366     return 0.
367     elif root.tagName == 'Altitude':
368     if root.hasAttribute("unit"):
369     unit=root.getAttribute("unit")
370     else:
371     unit="M"
372     for node in root.childNodes:
373     if isinstance(node, minidom.Text):
374 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
375 gross 977 return 0.
376     elif root.tagName == 'LongitudeLength':
377     if root.hasAttribute("unit"):
378     unit=root.getAttribute("unit")
379     else:
380     unit="M"
381     for node in root.childNodes:
382     if isinstance(node, minidom.Text):
383 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
384 gross 977 return 0.
385     elif root.tagName == 'LatitudeLength':
386     if root.hasAttribute("unit"):
387     unit=root.getAttribute("unit")
388     else:
389     unit="M"
390     for node in root.childNodes:
391     if isinstance(node, minidom.Text):
392 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
393 gross 977 return 0.
394     elif root.tagName == 'Depth':
395     if root.hasAttribute("unit"):
396     unit=root.getAttribute("unit")
397     else:
398     unit="M"
399     for node in root.childNodes:
400     if isinstance(node, minidom.Text):
401 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
402 gross 977 return 0.
403     elif root.tagName == 'LongitudeLength':
404     if root.hasAttribute("unit"):
405     unit=root.getAttribute("unit")
406     else:
407     unit="M"
408     for node in root.childNodes:
409     if isinstance(node, minidom.Text):
410 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
411 gross 977 return 0.
412     elif root.tagName == 'Tilt':
413     if root.hasAttribute("unit"):
414     unit=root.getAttribute("unit")
415     else:
416     unit="D.MM"
417     for node in root.childNodes:
418     if isinstance(node, minidom.Text):
419 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
420 gross 977 return 0.
421     elif root.tagName == 'Name':
422     for node in root.childNodes:
423     if isinstance(node, minidom.Text):
424     return node.nodeValue.strip()
425     return "unknown"
426 gross 981 elif root.tagName == 'Year':
427     for node in root.childNodes:
428     if isinstance(node, minidom.Text):
429     return int(node.nodeValue.strip())
430     return 0
431     elif root.tagName == 'Material':
432     for node in root.childNodes:
433     if isinstance(node, minidom.Text):
434     return node.nodeValue.strip()
435     return "unknown"
436     elif root.tagName == 'Extraction':
437     if root.hasAttribute("unit"):
438     unit=root.getAttribute("unit")
439     else:
440     unit="MT"
441     for node in root.childNodes:
442     if isinstance(node, minidom.Text):
443     return unitConverter(node.nodeValue.strip(), unit, "kg")
444     return 0.
445 gross 977 for node in root.childNodes:
446     if isinstance(node, minidom.Element): return _parse(node)
447    
448    
449     def parse(xml):
450     if isinstance(xml,str):
451     dom=minidom.parseString(xml)
452     else:
453     dom=minidom.parse(xml)
454     root=dom.getElementsByTagName('ESys')[0]
455     return _parse(root)
456    
457    
458     if __name__ == "__main__":
459    
460     FILE="newcastle_mining.xml"
461     mine=parse(open(FILE,'r'))
462     dsgn=Design(element_size=mine.getDiameter()*.2)
463     mine.fillDesign(dsgn)
464     dsgn.setScriptFileName("newcastle_mines.geo")
465     dsgn.setMeshFileName("newcastle_mines.msh")
466     print dsgn.getCommandString()
467 gross 981 print "start of records = ",mine.getStartOfRecords()
468 gross 977 dsgn.getMeshHandler()

  ViewVC Help
Powered by ViewVC 1.1.26