root/MGET/Branches/Jason/PythonPackage/src/GeoEco/Connectivity/CoralReefConnectivity.py @ 788

Revision 788, 87.7 KB (checked in by jjr8, 2 years ago)

Continued work on coral reef connectivity tool.
Line 
1# CoralReefConnectivity.py - Implements Eric Treml's coral reef
2# connectivity analysis.
3#
4# Copyright (C) 2008 Jason J. Roberts and Eric A. Treml
5#
6# This program is free software; you can redistribute it and/or
7# modify it under the terms of the GNU General Public License
8# as published by the Free Software Foundation; either version 2
9# of the License, or (at your option) any later version.
10#
11# This program is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14# GNU General Public License (available in the file LICENSE.TXT)
15# for more details.
16#
17# You should have received a copy of the GNU General Public License
18# along with this program; if not, write to the Free Software
19# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
20
21import os
22
23from GeoEco.DataProducts.Aviso import AvisoGriddedGeostrophicCurrents, _AvisoGriddedProduct_LongDescription
24from GeoEco.DataProducts.HYCOM import HYCOMGOMl0044D, _HYCOMGOMl004_LongDescription
25from GeoEco.DynamicDocString import DynamicDocString
26from GeoEco.Internationalization import _
27from GeoEco.Logging import Logger, ProgressReporter
28from GeoEco.Types import EnvelopeTypeMetadata
29
30class CoralReefConnectivity(object):
31    __doc__ = DynamicDocString()
32
33    @classmethod
34    def CreateSimulationFromArcGISRasters(cls, simulationDirectory, reefIDsRaster, reefCoverRaster, waterMaskRaster, crosses180=False, overwriteExisting=False):
35        cls.__doc__.Obj.ValidateMethodInvocation()
36
37        # Validate the coordinate systems, extents, and cell sizes of
38        # the input rasters.
39
40        from GeoEco.ArcGIS import GeoprocessorManager
41        gp = GeoprocessorManager.GetWrappedGeoprocessor()
42
43        describeReefIDsRaster = gp.Describe(reefIDsRaster)
44        describereefCoverRaster = gp.Describe(reefCoverRaster)
45        describeWaterMaskRaster = gp.Describe(waterMaskRaster)
46
47        if describeReefIDsRaster.SpatialReference.Type.lower() != u'projected':
48            Logger.RaiseException(ValueError(_(u'ArcGIS reports that the reef IDs raster %(raster)s uses a %(type)s coordinate system. It must use a projected coordinate system. Please project it and try again.') % {u'raster': reefIDsRaster, u'type': describeReefIDsRaster.SpatialReference.Type.lower()}))
49        if describeReefIDsRaster.SpatialReference.LinearUnitName.lower() != u'meter':
50            Logger.RaiseException(ValueError(_(u'ArcGIS reports that the reef IDs raster %(raster)s uses the linear unit "%(unit)s". It must use a coordinate system that has meters as its linear unit.  Please project it to a coordinate system that uses meters and try again.') % {u'raster': reefIDsRaster, u'unit': describeReefIDsRaster.SpatialReference.LinearUnitName.lower()}))
51
52        reefIDsRasterCS = gp.CreateSpatialReference(describeReefIDsRaster.SpatialReference).split(';')[0]
53        reefCoverRasterCS = gp.CreateSpatialReference(describereefCoverRaster.SpatialReference).split(';')[0]
54        waterMaskRasterCS = gp.CreateSpatialReference(describeWaterMaskRaster.SpatialReference).split(';')[0]
55
56        if reefCoverRasterCS.lower() != reefIDsRasterCS.lower():
57            Logger.RaiseException(ValueError(_(u'The reef cover raster %(raster1)s uses a different coordinate system than the reef IDs raster, %(raster2)s. Please project the reef cover raster to the reef IDs raster\'s coordinate system and try again.') % {u'raster1': reefCoverRaster, u'raster2': reefIDsRaster}))
58        if waterMaskRasterCS.lower() != reefIDsRasterCS.lower():
59            Logger.RaiseException(ValueError(_(u'The water mask raster %(raster1)s uses a different coordinate system than the reef IDs raster, %(raster2)s. Please project the water mask raster to the reef IDs raster\'s coordinate system and try again.') % {u'raster1': waterMaskRaster, u'raster2': reefIDsRaster}))
60
61        reefIDsRasterLeft, reefIDsRasterBottom, reefIDsRasterRight, reefIDsRasterTop = EnvelopeTypeMetadata.ParseFromArcGISString(describeReefIDsRaster.Extent)
62        reefCoverRasterLeft, reefCoverRasterBottom, reefCoverRasterRight, reefCoverRasterTop = EnvelopeTypeMetadata.ParseFromArcGISString(describereefCoverRaster.Extent)
63        waterMaskRasterLeft, waterMaskRasterBottom, waterMaskRasterRight, waterMaskRasterTop = EnvelopeTypeMetadata.ParseFromArcGISString(describeWaterMaskRaster.Extent)
64
65        if abs(reefCoverRasterLeft - reefIDsRasterLeft) > 0.001 or abs(reefCoverRasterBottom - reefIDsRasterBottom) > 0.001 or abs(reefCoverRasterRight - reefIDsRasterRight) > 0.001 or abs(reefCoverRasterTop - reefIDsRasterTop) > 0.001 or abs(describereefCoverRaster.MeanCellWidth - describeReefIDsRaster.MeanCellWidth) > 0.001 or abs(describereefCoverRaster.MeanCellHeight - describeReefIDsRaster.MeanCellHeight) > 0.001:
66            Logger.RaiseException(ValueError(_(u'The reef cover raster %(raster1)s has a different extent or cell size than the reef IDs raster, %(raster2)s. Please prepare a reef cover raster that has the same extent and cell size as the reef IDs raster and try again.') % {u'raster1': reefCoverRaster, u'raster2': reefIDsRaster}))
67        if abs(waterMaskRasterLeft - reefIDsRasterLeft) > 0.001 or abs(waterMaskRasterBottom - reefIDsRasterBottom) > 0.001 or abs(waterMaskRasterRight - reefIDsRasterRight) > 0.001 or abs(waterMaskRasterTop - reefIDsRasterTop) > 0.001 or abs(describeWaterMaskRaster.MeanCellWidth - describeReefIDsRaster.MeanCellWidth) > 0.001 or abs(describeWaterMaskRaster.MeanCellHeight - describeReefIDsRaster.MeanCellHeight) > 0.001:
68            Logger.RaiseException(ValueError(_(u'The water mask raster %(raster1)s has a different extent or cell size than the reef IDs raster, %(raster2)s. Please prepare a water mask raster that has the same extent and cell size as the reef IDs raster and try again.') % {u'raster1': reefCoverRaster, u'raster2': waterMaskRaster}))
69
70        # Validate that 0 is not used as a reef ID. The MATLAB code
71        # uses 0 to represent cells where no reefs are present.
72
73        from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
74        from GeoEco.DataManagement.Directories import Directory, TemporaryDirectory
75        import numpy
76
77        tempDir = TemporaryDirectory()
78        reefIDsImage, reefIDsNoDataValue = ArcGISRaster.ToNumpyArray(reefIDsRaster, tempRasterPath=os.path.join(tempDir.Path, u'reef_ids'))
79        if numpy.any(reefIDsImage == 0) and reefIDsNoDataValue != 0:
80            Logger.RaiseException(ValueError(_(u'The reef IDs raster %(raster)s uses 0 as a reef ID. This is not allowed.  Please remove the value 0 from the reef IDs raster and try again.')  % {u'raster': reefIDsRaster}))
81
82        # Validate that reef cover raster ranges between 0.0 and 1.0.
83
84        reefCoverImage, noDataValue = ArcGISRaster.ToNumpyArray(reefCoverRaster, tempRasterPath=os.path.join(tempDir.Path, u'reef_areas'))
85        if noDataValue is not None:
86            reefCoverImage[reefCoverImage == noDataValue] = 0
87        if numpy.any(reefCoverImage < 0) or numpy.any(reefCoverImage > 1):
88            Logger.RaiseException(ValueError(_(u'The reef cover raster %(raster)s includes values that are less than 0 or are greater than 1. The values of a cell of this raster is supposed to represent the proportion of the cell\'s area that is occupied by reef, thus the value is supposed to be between 0 and 1. Please prepare a raster that has the correct values and try again.')  % {u'raster': reefCoverRaster}))
89
90        # Create the simulation directory and the ReefData
91        # subdirectory.
92
93        oldLogInfoAsDebug = Logger.LogInfoAndSetInfoToDebug(_(u'Creating and initializing the simulation directory %(dir)s...') % {u'dir': simulationDirectory})
94        Logger.SetLogInfoAsDebug(True)
95        try:
96            from GeoEco.DataManagement.Directories import Directory, TemporaryDirectory
97
98            Directory.Create(simulationDirectory)
99            Directory.Create(os.path.join(simulationDirectory, u'ReefData'))
100
101            # Copy the reef IDs raster into the ReefData directory.
102
103            ArcGISRaster.Copy(reefIDsRaster, os.path.join(simulationDirectory, u'ReefData', u'reef_ids'))
104
105            # Copy the reef cover raster into the ReefData directory,
106            # setting the NoData values and non-reef cells to 0 in the process.
107
108            if GeoprocessorManager.GetArcGISMajorVersion() >= 10:
109                mapAlgebraExpression = u'float(con(isnull( [%(raster1)s] ) == 1 || isnull( [%(raster2)s] ) == 1, 0.0, [%(raster1)s] ))' % {u'raster1': reefCoverRaster, u'raster2': reefIDsRaster}
110            else:
111                mapAlgebraExpression = u'float(con(isnull( %(raster1)s ) == 1 || isnull( %(raster2)s ) == 1, 0.0, %(raster1)s ))' % {u'raster1': reefCoverRaster, u'raster2': reefIDsRaster}
112            gp.SingleOutputMapAlgebra_sa(mapAlgebraExpression, os.path.join(tempDir.Path, u'reef_areas2'))
113            ArcGISRaster.Copy(os.path.join(tempDir.Path, u'reef_areas2'), os.path.join(simulationDirectory, u'ReefData', u'reef_areas'))
114
115            # Copy the water mask raster into the ReefData directory,
116            # normalizing it to integer values where 1 is water, 0 is
117            # land.
118
119            if GeoprocessorManager.GetArcGISMajorVersion() >= 10:
120                mapAlgebraExpression = u'con(isnull( [%(raster)s] ) || [%(raster)s] == 0, 0, 1)' % {u'raster': waterMaskRaster}
121            else:
122                mapAlgebraExpression = u'con(isnull( %(raster)s ) || %(raster)s == 0, 0, 1)' % {u'raster': waterMaskRaster}
123            gp.SingleOutputMapAlgebra_sa(mapAlgebraExpression, os.path.join(tempDir.Path, u'water_mask'))
124            ArcGISRaster.Copy(os.path.join(tempDir.Path, u'water_mask'), os.path.join(simulationDirectory, u'ReefData', u'water_mask'))
125
126            # Create the reef_geometry file by calling Spatial
127            # Analyst's Zonal Geometry As Table tool. (Note that this
128            # tool will create a .dbf file no matter what, even if we
129            # give it a .csv or .txt extension. This is very annoying,
130            # because we can't read it easily without going through a
131            # database API, which is slow.)
132
133            gp.ZonalGeometryAsTable_sa(os.path.join(simulationDirectory, u'ReefData', u'reef_ids'), u'Value', os.path.join(simulationDirectory, u'ReefData', u'reef_geometry.dbf'), describeReefIDsRaster.MeanCellWidth)
134
135            # Create the directories that will hold the currents
136            # rasters.
137
138            Directory.Create(os.path.join(simulationDirectory, u'Currents'))
139            Directory.Create(os.path.join(simulationDirectory, u'Currents', u'u'))
140            Directory.Create(os.path.join(simulationDirectory, u'Currents', u'v'))
141
142            # Create the config file that stores properties of the
143            # simulation.
144
145            from ConfigParser import SafeConfigParser
146
147            scp = SafeConfigParser()
148            scp.add_section(u'Simulation')
149            scp.set(u'Simulation', u'Crosses180', unicode(crosses180))
150            scp.set(u'Simulation', u'CurrentsLoaded', unicode(False))
151            scp.set(u'Simulation', u'CurrentsProduct', u'')
152            f = file(os.path.join(simulationDirectory, u'Simulation.ini'), 'w')
153            try:
154                scp.write(f)
155            finally:
156                try:
157                    f.close()
158                except:
159                    pass
160
161        finally:
162            Logger.SetLogInfoAsDebug(oldLogInfoAsDebug)
163
164    @classmethod
165    def LoadAvisoGeostrophicCurrentsIntoSimulation(cls, simulationDirectory, username, password, product, startDate, endDate, timeout=60, maxRetryTime=120, cacheDirectory=None):
166        cls.__doc__.Obj.ValidateMethodInvocation()
167
168        # Parse and validate the Simulation.ini file.
169
170        scp, crosses180, currentsLoaded, currentsProduct = cls._ReadCurrentsInfoFromSimulationINI(simulationDirectory)
171
172        # If the simulation already has currents loaded in it from a
173        # different product, report an error.
174
175        if currentsLoaded and currentsProduct != u'Aviso ' + product:
176            Logger.RaiseException(ValueError(_(u'Cannot load Aviso %(prod1)s currents data into the simulation in directory %(dir)s because that simulation already has %(prod2)s currents data loaded into it. You can load additional %(prod2)s data, if you like, but to use Aviso %(prod1)s data, you must create a new simulation.') % {u'dir': simulationDirectory, u'prod1': product, u'prod2': currentsProduct}))
177
178        # First, download and create Aviso current rasters in an
179        # Aviso-specific directory. These will have the Aviso
180        # projection, extent, and cell size. We still need to project
181        # them to the simulation's projection, extent, and cell size.
182        #
183        # TODO: Need to handle rotation?
184
185        avisoCurrentsDirectory = os.path.join(simulationDirectory, 'OriginalAvisoCurrents')
186        try:
187            AvisoGriddedGeostrophicCurrents.CreateArcGISRasters(username, password, product, u'u', avisoCurrentsDirectory, rasterNameExpressions=['u', '%%Y', 'u%%Y%%j0000.img'], startDate=startDate, endDate=endDate, timeout=timeout, maxRetryTime=maxRetryTime, cacheDirectory=cacheDirectory)
188            AvisoGriddedGeostrophicCurrents.CreateArcGISRasters(username, password, product, u'v', avisoCurrentsDirectory, rasterNameExpressions=['v', '%%Y', 'v%%Y%%j0000.img'], startDate=startDate, endDate=endDate, timeout=timeout, maxRetryTime=maxRetryTime, cacheDirectory=cacheDirectory)
189
190        # If we created any rasters, record in the Simulation.ini file
191        # that we got some, so that any new rasters loaded into the
192        # simulation must be from the same product and depth.
193           
194        finally:
195            import glob
196            if len(glob.glob(os.path.join(os.path.join(avisoCurrentsDirectory, u'u', u'*', u'u[0-9][0-9][0-9][0-9][0-9][0-9][0-9]0000.img')))) > 0:
197                scp.set(u'Simulation', u'CurrentsLoaded', unicode(True))
198                scp.set(u'Simulation', u'CurrentsProduct', u'Aviso ' + product)
199                scp.set(u'Simulation', u'CurrentsDateType', u'Center')
200                if u'DT-Ref' in product or u'DT-Upd' in product:
201                    scp.set(u'Simulation', u'MaxSecondsBetweenCurrentsImages', unicode(7*86400))
202                else:
203                    scp.set(u'Simulation', u'MaxSecondsBetweenCurrentsImages', unicode(86400))
204                f = file(os.path.join(simulationDirectory, u'Simulation.ini'), 'w')
205                try:
206                    scp.write(f)
207                finally:
208                    try:
209                        f.close()
210                    except:
211                        pass
212
213        # Obtain the projection, extent, and cell size of the reef IDs
214        # raster, so we can project, snap, and clip the currents
215        # rasters to match it.
216
217        from GeoEco.ArcGIS import GeoprocessorManager
218        gp = GeoprocessorManager.GetWrappedGeoprocessor()
219
220        describeReefIDsRaster = gp.Describe(os.path.join(simulationDirectory, u'ReefData', u'reef_ids'))
221        reefIDsRasterCS = gp.CreateSpatialReference(describeReefIDsRaster.SpatialReference).split(';')[0]
222        reefIDsRasterLeft, reefIDsRasterBottom, reefIDsRasterRight, reefIDsRasterTop = EnvelopeTypeMetadata.ParseFromArcGISString(describeReefIDsRaster.Extent)
223        reefIDsCellSize = describeReefIDsRaster.MeanCellWidth
224
225        # Use a GeoEco function that calls Project_management to
226        # project, snap, and clip the rasters all in one step.
227
228        from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
229
230        oldGPExtent = gp.Extent
231        gp.Extent = describeReefIDsRaster.Extent
232        try:
233            ArcGISRaster.FindAndProjectClipAndOrExecuteMapAlgebra(avisoCurrentsDirectory,
234                                                                  os.path.join(simulationDirectory, u'Currents'),
235                                                                  projectedCoordinateSystem=reefIDsRasterCS,
236                                                                  resamplingTechnique=u'BILINEAR',
237                                                                  projectedCellSize=reefIDsCellSize,
238                                                                  registrationPoint=u'%g %g' % (reefIDsRasterLeft, reefIDsRasterBottom),
239                                                                  mapAlgebraExpression=u'setnull( inputRaster > 1000000 || inputRaster < -1000000, inputRaster / 100)',       # Convert from cm/s to to m/s and mask bogus values introduced by ArcGIS
240                                                                  wildcard=u'*.img',
241                                                                  searchTree=True,
242                                                                  outputRasterPythonExpression=u'os.path.join(outputWorkspace, os.path.splitext(inputRaster[len(workspaceToSearch)+1:])[0])',
243                                                                  modulesToImport=[u'os'],
244                                                                  skipExisting=True)
245
246        # Reset the geoprocessing extent environment setting to
247        # whatever the caller had before.
248
249        finally:
250            if oldGPExtent is not None:
251                gp.Extent = oldGPExtent
252
253        # Return successfully.
254
255        return simulationDirectory
256
257    @classmethod
258    def LoadHYCOMGOMl0044DCurrentsIntoSimulation(cls, simulationDirectory, startDate, endDate, depth=0., timeout=600, maxRetryTime=None, cacheDirectory=None):
259        cls.__doc__.Obj.ValidateMethodInvocation()
260
261        # Parse and validate the Simulation.ini file.
262
263        scp, crosses180, currentsLoaded, currentsProduct = cls._ReadCurrentsInfoFromSimulationINI(simulationDirectory)
264
265        # If the simulation already has currents loaded in it from a
266        # different product or depth, report an error.
267
268        if currentsLoaded and currentsProduct != u'HYCOM GOMl0.04 ' + unicode(int(depth)) + u' m':
269            Logger.RaiseException(ValueError(_(u'Cannot load %(prod1)s currents data into the simulation in directory %(dir)s because that simulation already has %(prod2)s currents data loaded into it. You can load additional %(prod2)s data, if you like, but to use %(prod1)s data, you must create a new simulation.') % {u'dir': simulationDirectory, u'prod1': u'HYCOM GOMl0.04 ' + unicode(int(depth)) + u' m', u'prod2': currentsProduct}))
270
271        # First, download and create HYCOM current rasters in a
272        # HYCOM-specific directory. These will have the HYCOM
273        # projection, extent, and cell size. We still need to project
274        # them to the simulation's projection, extent, and cell size.
275
276        hycomCurrentsDirectory = os.path.join(simulationDirectory, 'OriginalHYCOMCurrents')
277        try:
278            HYCOMGOMl0044D.CreateArcGISRasters(u'u', hycomCurrentsDirectory, rasterNameExpressions=['%(VariableName)s', '%%Y', '%(VariableName)s%%Y%%j0000.img'], minDepth=depth, maxDepth=depth, startDate=startDate, endDate=endDate, timeout=timeout, maxRetryTime=maxRetryTime, cacheDirectory=cacheDirectory)
279            HYCOMGOMl0044D.CreateArcGISRasters(u'v', hycomCurrentsDirectory, rasterNameExpressions=['%(VariableName)s', '%%Y', '%(VariableName)s%%Y%%j0000.img'], minDepth=depth, maxDepth=depth, startDate=startDate, endDate=endDate, timeout=timeout, maxRetryTime=maxRetryTime, cacheDirectory=cacheDirectory)
280
281        # If we created any rasters, record in the Simulation.ini file
282        # that we got some, so that any new rasters loaded into the
283        # simulation must be from the same product and depth.
284           
285        finally:
286            import glob
287            if len(glob.glob(os.path.join(os.path.join(hycomCurrentsDirectory, u'u', u'*', u'u[0-9][0-9][0-9][0-9][0-9][0-9][0-9]0000.img')))) > 0:
288                scp.set(u'Simulation', u'CurrentsLoaded', unicode(True))
289                scp.set(u'Simulation', u'CurrentsProduct', u'HYCOM GOMl0.04 ' + unicode(int(depth)) + u' m')
290                scp.set(u'Simulation', u'CurrentsDateType', u'Center')
291                scp.set(u'Simulation', u'MaxSecondsBetweenCurrentsImages', unicode(86400))
292                f = file(os.path.join(simulationDirectory, u'Simulation.ini'), 'w')
293                try:
294                    scp.write(f)
295                finally:
296                    try:
297                        f.close()
298                    except:
299                        pass
300
301        # Obtain the projection, extent, and cell size of the reef IDs
302        # raster, so we can project, snap, and clip the currents
303        # rasters to match it.
304
305        from GeoEco.ArcGIS import GeoprocessorManager
306        gp = GeoprocessorManager.GetWrappedGeoprocessor()
307
308        describeReefIDsRaster = gp.Describe(os.path.join(simulationDirectory, u'ReefData', u'reef_ids'))
309        reefIDsRasterCS = gp.CreateSpatialReference(describeReefIDsRaster.SpatialReference).split(';')[0]
310        reefIDsRasterLeft, reefIDsRasterBottom, reefIDsRasterRight, reefIDsRasterTop = EnvelopeTypeMetadata.ParseFromArcGISString(describeReefIDsRaster.Extent)
311        reefIDsCellSize = describeReefIDsRaster.MeanCellWidth
312
313        # Use a GeoEco function that calls Project_management to
314        # project, snap, and clip the rasters all in one step.
315
316        from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
317
318        oldGPExtent = gp.Extent
319        gp.Extent = describeReefIDsRaster.Extent
320        try:
321            ArcGISRaster.FindAndProjectClipAndOrExecuteMapAlgebra(hycomCurrentsDirectory,
322                                                                  os.path.join(simulationDirectory, u'Currents'),
323                                                                  projectedCoordinateSystem=reefIDsRasterCS,
324                                                                  resamplingTechnique=u'BILINEAR',
325                                                                  projectedCellSize=reefIDsCellSize,
326                                                                  registrationPoint=u'%g %g' % (reefIDsRasterLeft, reefIDsRasterBottom),
327                                                                  wildcard=u'*.img',
328                                                                  searchTree=True,
329                                                                  outputRasterPythonExpression=u'os.path.join(outputWorkspace, os.path.splitext(inputRaster[len(workspaceToSearch)+1:])[0])',
330                                                                  modulesToImport=[u'os'],
331                                                                  skipExisting=True)
332
333        # Reset the geoprocessing extent environment setting to
334        # whatever the caller had before.
335
336        finally:
337            if oldGPExtent is not None:
338                gp.Extent = oldGPExtent
339
340        # Return successfully.
341
342        return simulationDirectory
343
344    @classmethod
345    def _ReadCurrentsInfoFromSimulationINI(cls, simulationDirectory):
346        if not os.path.isfile(os.path.join(simulationDirectory, u'Simulation.ini')):
347            Logger.RaiseException(ValueError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: it does not contain a file called Simulation.ini. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory}))
348
349        from ConfigParser import SafeConfigParser
350
351        scp = SafeConfigParser()
352        f = file(os.path.join(simulationDirectory, u'Simulation.ini'), 'r')
353        try:
354            try:
355                scp.readfp(f, os.path.join(simulationDirectory, u'Simulation.ini'))
356            except:
357                Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: the file Simulation.ini in that directory could not be parsed. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
358                raise
359        finally:
360            try:
361                f.close()
362            except:
363                pass
364
365        try:
366            crosses180 = scp.getboolean(u'Simulation', u'Crosses180')
367        except:
368            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a boolean option named Crosses180 from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
369            raise
370
371        try:
372            currentsLoaded = scp.getboolean(u'Simulation', u'CurrentsLoaded')
373        except:
374            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a boolean option named CurrentsLoaded from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
375            raise
376
377        try:
378            currentsProduct = unicode(scp.get(u'Simulation', u'CurrentsProduct'))
379        except:
380            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a string option named CurrentsProduct from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
381            raise
382
383        return scp, crosses180, currentsLoaded, currentsProduct
384
385    @classmethod
386    def RunSimulation(cls, simulationDirectory, outputDirectory, startDate, duration=30.0, simulationTimeStep=2.4, summarizationPeriod=10, initialLarvaeDensity=10000.0, densityRasterCutoff=0.1, diffusivity=25.0, includeReefIDs=None, excludeReefIDs=None, overwriteExisting=False):
387        cls.__doc__.Obj.ValidateMethodInvocation()
388
389        # Parse and validate the Simulation.ini file.
390
391        if not os.path.isfile(os.path.join(simulationDirectory, u'Simulation.ini')):
392            Logger.RaiseException(ValueError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: it does not contain a file called Simulation.ini. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory}))
393
394        from ConfigParser import SafeConfigParser
395
396        scp = SafeConfigParser()
397        f = file(os.path.join(simulationDirectory, u'Simulation.ini'), 'r')
398        try:
399            try:
400                scp.readfp(f, os.path.join(simulationDirectory, u'Simulation.ini'))
401            except:
402                Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: the file Simulation.ini in that directory could not be parsed. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
403                raise
404        finally:
405            try:
406                f.close()
407            except:
408                pass
409
410        try:
411            crosses180 = scp.getboolean(u'Simulation', u'Crosses180')
412        except:
413            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a boolean option named Crosses180 from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
414            raise
415
416        try:
417            currentsLoaded = scp.getboolean(u'Simulation', u'CurrentsLoaded')
418        except:
419            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a boolean option named CurrentsLoaded from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool and try again.') % {u'dir': simulationDirectory})
420            raise
421        if not currentsLoaded:
422            Logger.RaiseException(ValueError(_(u'The coral reef connectivity simulation in directory %(dir)s does not contain any ocean currents data. Please load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory}))
423
424        try:
425            currentsProduct = unicode(scp.get(u'Simulation', u'CurrentsProduct'))
426        except:
427            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a string option named CurrentsProduct from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool, load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory})
428            raise
429
430        try:
431            currentsDateType = unicode(scp.get(u'Simulation', u'CurrentsDateType'))
432        except:
433            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse a string option named CurrentsProduct from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool, load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory})
434            raise
435        if currentsDateType.lower() not in [u'center']:
436            Logger.RaiseException(ValueError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: the CurrentsDataType option in the Simulation.ini file has the unknown value %(val)s. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool, load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory, u'val': currentsDateType}))
437
438        try:
439            maxSecondsBetweenCurrentsImages = scp.getint(u'Simulation', u'MaxSecondsBetweenCurrentsImages')
440        except:
441            Logger.LogExceptionAsError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: failed to parse an integer option named MaxSecondsBetweenCurrentsImages from the file Simulation.ini in that directory. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool, load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory})
442            raise
443        if maxSecondsBetweenCurrentsImages <= 0:
444            Logger.RaiseException(ValueError(_(u'The directory %(dir)s does not appear to be a properly-initialized coral reef connectivity simulation directory: the MaxSecondsBetweenCurrentsImages option in the Simulation.ini file is less than or equal to zero. It must be greater than zero. Please create a simulation directory using the Create Coral Reef Connectivity Simulation tool, load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory}))
445
446        # Validate other input parameters.
447
448        if duration <= 0:
449            Logger.RaiseException(ValueError(_(u'The simulation duration must be greater than zero.')))
450
451        if duration - simulationTimeStep/24 <= 0:
452            Logger.RaiseException(ValueError(_(u'The time step must be shorter than or equal to the simulation duration. For accurate results, the time step should be significantly shorter than the simulation duration. For example, if the simulation duration is 60 days, a time step of 1 hour would be appropriate.')))
453
454        if includeReefIDs is not None and excludeReefIDs is not None:
455            Logger.RaiseException(ValueError(_(u'You cannot specify both a list of reefs to include and a list of reefs to exclude. You must specify one, or the other, or neither.')))
456
457        # Build lists of the currents rasters that are loaded into the
458        # simulation, and validate that we have currents for the start
459        # date and duration specified by the caller.
460
461        import datetime
462        import glob
463
464        uRasters = glob.glob(os.path.join(simulationDirectory, u'Currents', u'u', u'*', u'u[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]'))
465        uRasters.sort()
466        vRasters = glob.glob(os.path.join(simulationDirectory, u'Currents', u'v', u'*', u'v[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]'))
467        vRasters.sort()
468
469        if len(uRasters) <= 0:
470            Logger.RaiseException(ValueError(_(u'The coral reef connectivity simulation in directory %(dir)s does not contain any ocean currents data. Please load ocean currents into it using a tool designed for this purpose, and try again.') % {u'dir': simulationDirectory}))
471        if len(uRasters) != len(vRasters):
472            Logger.RaiseException(ValueError(_(u'The ocean currents data in the coral reef connectivity simulation in directory %(dir)s appears to be incompletely loaded. The number of "u" rasters does not equal the number of "v" rasters, indicating that the load operation did not complete successfully. Please try loading ocean currents again, and then try to run the simulation.') % {u'dir': simulationDirectory}))
473
474        rasterDates = []
475
476        for i in range(len(uRasters)):
477            uRasterName = os.path.basename(uRasters[i])
478            if uRasterName[1:] != os.path.basename(vRasters[i])[1:]:
479                Logger.RaiseException(ValueError(_(u'The ocean currents data in the coral reef connectivity simulation in directory %(dir)s appears to be incompletely loaded. The "u" raster %(r1)s could not be matched up with a "v" raster with the same date (the next available "v" raster is %(r2)s). Please try loading ocean currents again, and then try to run the simulation.') % {u'dir': simulationDirectory, u'r1': uRasters[i], u'r2': vRasters[i]}))
480            rasterDates.append(datetime.datetime(int(uRasterName[1:5]), 1, 1) + datetime.timedelta(days=int(uRasterName[5:8]) - 1, hours=int(uRasterName[8:10]), minutes=int(uRasterName[10:12])))
481
482        if currentsDateType.lower() == u'center':
483            currentsDateStartDelta = datetime.timedelta(seconds=maxSecondsBetweenCurrentsImages / 2)
484            currentsDateEndDelta = datetime.timedelta(seconds=maxSecondsBetweenCurrentsImages / 2)
485        else:
486            Logger.RaiseException(NotImplementedError(_(u'This tool does not currently support a CurrentsDateType of "%(type)s". Please contact the author of this tool for assistance.') % {u'type': currentsDateType}))
487
488        if startDate < rasterDates[0] - currentsDateStartDelta:
489            Logger.RaiseException(ValueError(_(u'The start date of the simulation (%(start)s) occurs too far before the date of the first ocean currents image (%(date)s) that is loaded in the coral reef connectivity simulation in directory %(dir)s. To fix this problem, either move the start date forward or load some older ocean currents data into the simulation, so that the start date matches up with the currents data.') % {u'dir': simulationDirectory, u'start': str(startDate), u'date': str(rasterDates[0])}))
490
491        if startDate > rasterDates[-1] + currentsDateEndDelta:
492            Logger.RaiseException(ValueError(_(u'The start date of the simulation (%(start)s) occurs too far after the date of the last ocean currents image (%(date)s) that is loaded in the coral reef connectivity simulation in directory %(dir)s. To fix this problem, either move the start date backward or load some more recent ocean currents data into the simulation, so that the start date matches up with the currents data.') % {u'dir': simulationDirectory, u'start': str(startDate), u'date': str(rasterDates[-1])}))
493
494        endDate = startDate + datetime.timedelta(days=duration)
495
496        if endDate > rasterDates[-1] + currentsDateEndDelta:
497            Logger.RaiseException(ValueError(_(u'The end date of the simulation (%(end)s) occurs too far after the date of the last ocean currents image (%(date)s) that is loaded in the coral reef connectivity simulation in directory %(dir)s. To fix this problem, either move the start date backward, reduce the duration of the simulation, or load some more recent ocean currents data into the simulation, so that the end date matches up with the currents data.') % {u'dir': simulationDirectory, u'end': str(endDate), u'date': str(rasterDates[-1])}))
498
499        startRasterIndex = 0
500        while startDate > rasterDates[startRasterIndex] + currentsDateEndDelta:
501            startRasterIndex += 1
502
503        endRasterIndex = startRasterIndex
504        while endDate > rasterDates[endRasterIndex] + currentsDateEndDelta:
505            endRasterIndex += 1
506            if rasterDates[endRasterIndex] - rasterDates[endRasterIndex-1] > datetime.timedelta(seconds=maxSecondsBetweenCurrentsImages):
507                Logger.RaiseException(ValueError(_(u'The ocean currents data that is loaded in the coral reef connectivity simulation in directory %(dir)s has a data gap in the range of dates between the simulation start date (%(start)s) and end date (%(end)s). A gap of %(gap)s occurs between %(d1)s and %(d2)s, which is larger than the maximum time permitted between images (%(max)s) for %(prod)s data. To fix this problem, either adjust the start date or duration, or load ocean currents data into the simulation that fills the gap.') % {u'dir': simulationDirectory, u'start': str(startDate), u'end': str(endDate), u'gap': str(rasterDates[endRasterIndex] - rasterDates[endRasterIndex-1]), u'd1': rasterDates[endRasterIndex-1], u'd2': rasterDates[endRasterIndex], u'max': datetime.timedelta(seconds=maxSecondsBetweenCurrentsImages), u'prod': currentsProduct}))
508
509        # Read the reef rasters and geometry table into 2D numpy
510        # arrays.
511
512        Logger.Info(_(u'Reading coral reef data...'))
513
514        from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
515
516        reefIDsImage, reefIDsNoDataValue = ArcGISRaster.ToNumpyArray(os.path.join(simulationDirectory, u'ReefData', u'reef_ids'))
517        reefIDsNoDataValue = int(reefIDsNoDataValue)
518        reefCoverImage = ArcGISRaster.ToNumpyArray(os.path.join(simulationDirectory, u'ReefData', u'reef_areas'))[0]
519        waterMaskImage = ArcGISRaster.ToNumpyArray(os.path.join(simulationDirectory, u'ReefData', u'water_mask'))[0]
520
521        from GeoEco.DatabaseAccess.ArcGIS import ArcGIS91DatabaseConnection
522        from GeoEco.DatabaseAccess.InMemory import InMemoryDatabaseConnection
523
524        arcGISConn = ArcGIS91DatabaseConnection()
525        inMemoryConn = InMemoryDatabaseConnection()
526       
527        if includeReefIDs is not None:
528            inMemoryConn.ImportTable(arcGISConn, os.path.join(simulationDirectory, u'ReefData', u'reef_geometry.dbf'), u'reef_geometry', [u'VALUE', u'XCENTROID', u'YCENTROID'], where=u'"VALUE" IN (%s)' % u', '.join(map(unicode, includeReefIDs)), orderBy=[u'VALUE'], directions=[u'Ascending'])
529            if inMemoryConn.GetRowCount(u'reef_geometry') <= 0:
530                Logger.RaiseException(ValueError(_(u'The list of reefs IDs to include in the simulation does contain any IDs that also exist in the simulation directory %(dir)s. Please specify at least one existing reef ID.') % {u'dir': simulationDirectory}))
531               
532        elif excludeReefIDs is not None:
533            inMemoryConn.ImportTable(arcGISConn, os.path.join(simulationDirectory, u'ReefData', u'reef_geometry.dbf'), u'reef_geometry', [u'VALUE', u'XCENTROID', u'YCENTROID'], where=u'"VALUE" NOT IN (%s)' % u', '.join(map(unicode, excludeReefIDs)), orderBy=[u'VALUE'], directions=[u'Ascending'])
534            if inMemoryConn.GetRowCount(u'reef_geometry') <= 0:
535                Logger.RaiseException(ValueError(_(u'The list of reefs IDs to exclude from the simulation excluded all of the reefs in the simulation directory %(dir)s. Please remove some IDs from this list so that at least one reef will be included in the simulation.') % {u'dir': simulationDirectory}))
536
537        else:
538            inMemoryConn.ImportTable(arcGISConn, os.path.join(simulationDirectory, u'ReefData', u'reef_geometry.dbf'), u'reef_geometry', [u'VALUE', u'XCENTROID', u'YCENTROID'], orderBy=[u'VALUE'], directions=[u'Ascending'])
539            if inMemoryConn.GetRowCount(u'reef_geometry') <= 0:
540                Logger.RaiseException(ValueError(_(u'The reef data in simulation directory %(dir)s does not contain any reefs. Please recreate the simulation using input rasters that contain at least one reef.') % {u'dir': simulationDirectory}))
541
542        # Read the ocean currents into parallel lists of 2D numpy
543        # arrays.
544
545        import numpy
546
547        imagesToRead = (endRasterIndex - startRasterIndex + 1) * 2
548        Logger.Info(_(u'Reading %i ocean currents images...') % imagesToRead)
549        progressReporter = ProgressReporter(progressMessage1=_(u'Still reading: %(elapsed)s elapsed, %(opsCompleted)i images read, %(perOp)s per image, %(opsRemaining)i remaining, estimated completion time: %(etc)s.'),
550                                            completionMessage=_(u'Finished reading: %(elapsed)s elapsed, %(opsCompleted)i images read, %(perOp)s per image.'))
551        progressReporter.Start(imagesToRead)
552
553        uImageList = []
554        vImageList = []
555        uvDateList = []
556
557        i = startRasterIndex
558        while i >= startRasterIndex and i <= endRasterIndex:
559            imageDate = rasterDates[i]
560            uvDateList.append(imageDate)
561           
562            image, noDataValue = ArcGISRaster.ToNumpyArray(os.path.join(simulationDirectory, u'Currents', u'u', unicode(imageDate.year), imageDate.strftime('u%Y%j%H%M')))
563            image[image == noDataValue] = numpy.nan
564            uImageList.append(image)
565            progressReporter.ReportProgress()
566           
567            image, noDataValue = ArcGISRaster.ToNumpyArray(os.path.join(simulationDirectory, u'Currents', u'v', unicode(imageDate.year), imageDate.strftime('v%Y%j%H%M')))
568            image[image == noDataValue] = numpy.nan
569            vImageList.append(image)
570            progressReporter.ReportProgress()
571           
572            i += 1
573
574        # Stack the numpy arrays into two 3D arrays that we will pass
575        # to the MATLAB function.
576        #
577        # Note that with numpy, it appears that 2D arrays are
578        # traditionally indexed [y,x] but that 3D arrays are [y,x,t].
579        # This is what is output by numpy's dstack function. This is
580        # kind of screwy, because when you print a 3D numpy array, it
581        # looks much better if ordered [t,y,x] than [y,x,t]. But
582        # MATLAB was the inspiration for numpy and [y,x,t] is
583        # traditional in MATLAB as well. Finally, Eric Treml's
584        # original MATLAB code used [y,x,z].
585
586        uImages = numpy.dstack(tuple(uImageList))
587        del uImageList
588       
589        vImages = numpy.dstack(tuple(vImageList))
590        del vImageList
591
592        # Build a list that specifies the t index into the 3D arrays
593        # for each time step. Note that when we pass this list to
594        # MATLAB, we must increment all of the indices by 1, because
595        # MATLAB uses 1-based indexing (Python uses 0-based).
596
597        import math       
598
599        numTimeSteps = int(math.ceil(duration / (simulationTimeStep/24)))
600        uvIndexForTimestep = [0]
601        for i in range(1, numTimeSteps):
602            if startDate + datetime.timedelta(hours=simulationTimeStep*i) <= uvDateList[uvIndexForTimestep[-1]] + currentsDateEndDelta:
603                uvIndexForTimestep.append(uvIndexForTimestep[-1])
604            else:
605                uvIndexForTimestep.append(uvIndexForTimestep[-1] + 1)
606
607        # Look up the cell size of the rasters.
608
609        from GeoEco.ArcGIS import GeoprocessorManager
610        gp = GeoprocessorManager.GetWrappedGeoprocessor()
611        describeReefIDsRaster = gp.Describe(os.path.join(simulationDirectory, u'ReefData', u'reef_ids'))
612        cellSize = describeReefIDsRaster.MeanCellWidth
613
614        # Calculate and report the maximum Courant number. Issue a
615        # warning if it is greater than or equal to 1.0 because the
616        # simulation is likely to be unstable. Provide an estimate of
617        # the largest time step that would allow the Courant number to
618        # be less than 1.0.
619
620        hasData = numpy.logical_and(numpy.logical_not(numpy.isnan(uImages)), numpy.logical_not(numpy.isnan(vImages)))
621        maxVelocity = max(numpy.max(uImages[hasData]), numpy.max(vImages[hasData]))
622        maxCourant = maxVelocity * (simulationTimeStep*3600) / cellSize
623
624        if maxCourant <= 0.25:
625            Logger.Info(_(u'The maximum Courant number is %(mc)f, which is less than or equal to 0.25. The simulation is likely to be numerically stable.') % {u'mc': maxCourant})
626        else:
627            maxTimeStep = cellSize / maxVelocity / 3600 * 0.25
628            if maxCourant <= 0.5:
629                Logger.Warning(_(u'The maximum Courant number is %(mc)f, which is greater than 0.25 and less than or equal to 0.5. The simulation may exhibit some instability. Please review the results carefully. To improve the chance that the simulation will be stable, we recommend you reduce the time step to %(mts)g or less, so that the maximum Courant number is less than or equal to 0.25.') % {u'mc': maxCourant, u'mts': maxTimeStep})
630            else:
631                Logger.Warning(_(u'The maximum Courant number is %(mc)f, which is greater than 0.5. The simulation is likely to be unstable. Please review the results carefully. To improve the chance that the simulation will be stable, we recommend you reduce the time step to %(mts)g or less, so that the maximum Courant number is less than or equal to 0.25.') % {u'mc': maxCourant, u'mts': maxTimeStep})
632
633        # Create a temporary directory and write the arrays to it in
634        # binary format.
635
636        reefIDs = inMemoryConn.GetFieldValues(u'reef_geometry', u'VALUE')
637
638        from GeoEco.DataManagement.Directories import TemporaryDirectory
639        tempDir = TemporaryDirectory()
640
641        try:
642            reefIDsArray = numpy.array(reefIDs)
643            reefIDsDataType = reefIDsArray.dtype.name
644            reefIDsFile = os.path.join(tempDir.Path, 'ReefIDs.bin')
645            reefIDsArray.tofile(reefIDsFile)
646
647            reefIDsImageDataType = reefIDsImage.dtype.name
648            reefIDsImageFile = os.path.join(tempDir.Path, 'ReefIDsImage.bin')
649            reefIDsImage.tofile(reefIDsImageFile)
650
651            reefCoverImageDataType = reefCoverImage.dtype.name
652            reefCoverImageFile = os.path.join(tempDir.Path, 'ReefCoverImage.bin')
653            reefCoverImage.tofile(reefCoverImageFile)
654
655            waterMaskImageDataType = waterMaskImage.dtype.name
656            waterMaskImageFile = os.path.join(tempDir.Path, 'waterMaskImage.bin')
657            waterMaskImage.tofile(waterMaskImageFile)
658
659            uImagesDataType = uImages.dtype.name
660            uImagesFile = os.path.join(tempDir.Path, 'uImages.bin')
661            uImages.tofile(uImagesFile)
662
663            vImagesDataType = vImages.dtype.name
664            vImagesFile = os.path.join(tempDir.Path, 'vImages.bin')
665            vImages.tofile(vImagesFile)
666
667            uvIndexForTimestepArray = numpy.array(uvIndexForTimestep) + 1
668            uvIndexForTimestepDataType = uvIndexForTimestepArray.dtype.name
669            uvIndexForTimestepFile = os.path.join(tempDir.Path, 'UVIndexForTimestep.bin')
670            uvIndexForTimestepArray.tofile(uvIndexForTimestepFile)
671
672            # Execute RunLarvalDispersal.py to run the simulation.
673            # This script calls MATLAB functions. We prefer to call
674            # those functions directly right here but there is a
675            # continuing incompatibility between MATLAB DLLs and
676            # ArcGIS DLLs (they both try to load their own
677            # incompatible versions of xerces-c_2_7.dll) so we have to
678            # do it in a separate process.
679
680            dispersalMatrixFile = os.path.join(tempDir.Path, u'DispersalMatrix.bin')
681            densityImagesFile = os.path.join(tempDir.Path, u'DensityImages.bin')
682
683            from GeoEco.DataManagement.Processes import ChildProcess
684            import win32api
685
686            ChildProcess.ExecuteProgram(win32api.FindExecutable(os.path.join(os.path.dirname(__file__), 'RunLarvalDispersal.py'), os.path.dirname(__file__))[1],
687                                        arguments=[os.path.join(os.path.dirname(__file__), 'RunLarvalDispersal.py'),
688                                                   repr(uImages.shape[2]),      # t
689                                                   repr(uImages.shape[0]),      # y
690                                                   repr(uImages.shape[1]),      # x
691                                                   reefIDsFile,
692                                                   reefIDsDataType,
693                                                   reefIDsImageFile,
694                                                   reefIDsImageDataType,
695                                                   reefCoverImageFile,
696                                                   reefCoverImageDataType,
697                                                   waterMaskImageFile,
698                                                   waterMaskImageDataType,
699                                                   uImagesFile,
700                                                   uImagesDataType,
701                                                   vImagesFile,
702                                                   vImagesDataType,
703                                                   repr(cellSize),
704                                                   repr(simulationTimeStep * 3600.0),
705                                                   repr(initialLarvaeDensity / 1000000.0),
706                                                   repr(summarizationPeriod),
707                                                   repr(diffusivity),
708                                                   uvIndexForTimestepFile,
709                                                   uvIndexForTimestepDataType,
710                                                   dispersalMatrixFile,
711                                                   densityImagesFile],
712                                        stdoutLogLevel=u'Info',
713                                        windowState=u'invisible',
714                                        maxRunTime=None)
715
716            del uImages, vImages
717
718            # Read the output files.
719
720            dispersalMatrix = numpy.fromfile(dispersalMatrixFile, 'float32')
721            dispersalMatrix = dispersalMatrix.reshape(len(reefIDs), len(reefIDs), len(dispersalMatrix) / len(reefIDs) / len(reefIDs))
722
723            densityImages = numpy.fromfile(densityImagesFile, 'float32')
724            densityImages = densityImages.reshape(reefIDsImage.shape[0], reefIDsImage.shape[1], len(densityImages) / reefIDsImage.shape[0] / reefIDsImage.shape[1])
725
726        finally:
727            del tempDir
728
729        # The DisperseLarvae function returns densities in particles
730        # per square meter. Convert to particles per square km.
731
732        densityImages *= 1000000.0
733
734        # Mask cells that are land and that have a density that is
735        # below the threshold.
736
737        densityImages[waterMaskImage == 0, :] = 0
738        if densityRasterCutoff is not None:
739            densityImages[densityImages < initialLarvaeDensity * densityRasterCutoff / 100.0] = 0
740
741        # Create the output personal geodatabase in the output
742        # directory.
743
744        outputGDB = os.path.join(outputDirectory, u'ConnectivityGeodatabase.mdb')
745        if gp.Exists(outputGDB):
746            if not overwriteExisting:
747                Logger.RaiseException(ValueError(_(u'The output geodatabase %s already exists. Please delete it or specify that existing outputs should be overwritten and try again.') % outputGDB))
748            gp.Delete_management(outputGDB)
749
750        gp.CreatePersonalGDB_management(outputDirectory, u'ConnectivityGeodatabase.mdb')
751
752        # Create the edge list feature class.
753
754        gp.CreateFeatureClass_management(outputGDB, u'Edges', u'POLYLINE', None, u'DISABLED', u'DISABLED', gp.CreateSpatialReference(describeReefIDsRaster.SpatialReference).split(';')[0])
755        gp.AddField_management(os.path.join(outputGDB, u'Edges'), u'FromReefID', u'LONG')
756        gp.AddField_management(os.path.join(outputGDB, u'Edges'), u'ToReefID', u'LONG')
757        gp.AddField_management(os.path.join(outputGDB, u'Edges'), u'MaxDispersal', u'FLOAT')
758
759        # Populate the edge list feature class.
760
761        maxDispersal = dispersalMatrix[:,:,1:].max(2)
762        for fromReef in range(len(reefIDs)):
763            for toReef in range(len(reefIDs)):
764                if fromReef != toReef:
765                    maxDispersal[fromReef, toReef] /= dispersalMatrix[fromReef, fromReef, 0]
766                else:
767                    maxDispersal[fromReef, toReef] = 0
768
769        nonZeroEdges = sum(sum(maxDispersal > 0.0001))
770        xCentroids = inMemoryConn.GetFieldValues(u'reef_geometry', u'XCENTROID')
771        yCentroids = inMemoryConn.GetFieldValues(u'reef_geometry', u'YCENTROID')
772        shapeFieldName = gp.Describe(os.path.join(outputGDB, u'Edges')).ShapeFieldName
773
774        if nonZeroEdges > 0:
775            Logger.Info(_(u'Writing %i edges to the edge list in the output geodatabase...') % nonZeroEdges)
776            cur = arcGISConn.OpenInsertCursor(os.path.join(outputGDB, u'Edges'), rowCount=nonZeroEdges)
777            for fromReef in range(len(reefIDs)):
778                for toReef in range(len(reefIDs)):
779                    if maxDispersal[fromReef, toReef] > 0.0001:
780                        point = gp.CreateObject(u'Point')
781                        line = gp.CreateObject(u'Array')
782                        point.X = xCentroids[fromReef]
783                        point.Y = yCentroids[fromReef]
784                        line.Add(point)
785                        point.X = xCentroids[toReef]
786                        point.Y = yCentroids[toReef]
787                        line.Add(point)
788                        cur.SetValue(shapeFieldName, line)
789                        cur.SetValue(u'FromReefID', reefIDs[fromReef])
790                        cur.SetValue(u'ToReefID', reefIDs[toReef])
791                        cur.SetValue(u'MaxDispersal', float(maxDispersal[fromReef, toReef]))        # Must convert from numpy float to Python float
792                        cur.InsertRow()
793        else:
794            Logger.Warning(_(u'The edge list in the output geodatabase will be empty because none of the reefs are connected.'))
795
796        # Create the DensityRasters subdirectory, if it does not
797        # already exist.
798       
799        from GeoEco.DataManagement.Directories import Directory
800        from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
801             
802        if not Directory.Exists(os.path.join(outputDirectory, u'DensityRasters'))[0]:
803            Directory.Create(os.path.join(outputDirectory, u'DensityRasters'))
804
805        # If the directory already exists and the caller requested
806        # that we overwrite existing outputs, delete any existing
807        # density rasters.
808
809        elif overwriteExisting:
810            oldLogInfoAsDebug = Logger.GetLogInfoAsDebug()
811            Logger.SetLogInfoAsDebug(True)
812            try:
813                ArcGISRaster.FindAndDelete(os.path.join(outputDirectory, u'DensityRasters'), u'*')
814            finally:
815                Logger.SetLogInfoAsDebug(oldLogInfoAsDebug)
816
817        # Create the density rasters in the subdirectory.
818
819        coordinateSystem = gp.CreateSpatialReference(describeReefIDsRaster.SpatialReference).split(';')[0]
820
821        Logger.Info(_(u'Writing %i density rasters to the output directory...') % densityImages.shape[2])
822        progressReporter = ProgressReporter(progressMessage1=_(u'Still writing density rasters: %(elapsed)s elapsed, %(opsCompleted)i rasters written, %(perOp)s per raster, %(opsRemaining)i remaining, estimated completion time: %(etc)s.'),
823                                            completionMessage=_(u'Finished writing density rasters: %(elapsed)s elapsed, %(opsCompleted)i rasters written, %(perOp)s per raster.'))
824        progressReporter.Start(densityImages.shape[2])
825
826        for i in range(densityImages.shape[2]):
827            oldLogInfoAsDebug = Logger.GetLogInfoAsDebug()
828            Logger.SetLogInfoAsDebug(True)
829            try:
830                ArcGISRaster.FromNumpyArray(densityImages[:,:,i].copy(),    # Copy is currently needed because the current implementation of ArcGISRaster.FromNumpyArray requires the array to be in C order.
831                                            os.path.join(outputDirectory, u'DensityRasters', (startDate + datetime.timedelta(seconds=simulationTimeStep*3600.0*summarizationPeriod*i)).strftime('d%Y%j%H%M')),
832                                            EnvelopeTypeMetadata.ParseFromArcGISString(describeReefIDsRaster.Extent)[0],
833                                            EnvelopeTypeMetadata.ParseFromArcGISString(describeReefIDsRaster.Extent)[1],
834                                            cellSize,
835                                            nodataValue=0,
836                                            coordinateSystem=coordinateSystem,
837                                            buildPyramids=True,
838                                            overwriteExisting=overwriteExisting)
839            finally:
840                Logger.SetLogInfoAsDebug(oldLogInfoAsDebug)
841
842            progressReporter.ReportProgress()
843
844        # Create a raster catalog in the output geodatabase for
845        # the density rasters, so they can be easily animated. Add
846        # and populate StartTime and EndTime fields, which define
847        # the period for which the raster applies, as recommended
848        # by the ArcGIS documentation.
849
850        Logger.Info(_(u'Creating a raster catalog for the density rasters...'))
851
852        gp.CreateRasterCatalog_management(outputGDB, u'DensityRasters', coordinateSystem, coordinateSystem, None, None, None, None, u'Unmanaged')
853        gp.AddField_management(os.path.join(outputGDB, u'DensityRasters'), u'StartTime', u'DATE')
854        gp.AddField_management(os.path.join(outputGDB, u'DensityRasters'), u'EndTime', u'DATE')
855        gp.WorkspaceToRasterCatalog(os.path.join(outputDirectory, u'DensityRasters'), os.path.join(outputGDB, u'DensityRasters'))
856
857        # Populate the StartTime and EndTime fields of the raster
858        # catalog. Originally, I implemented this using the ArcGIS
859        # Calculate Fields tool, but it kept crashing in 9.3,
860        # regardless of whether I used Python or VB, and regardless of
861        # whether I created the geoprocessor using
862        # arcgisscripting.create or win32com.client.Dispatch. To work
863        # around this ArcGIS bug, I'm using the MGET Calculate Fields
864        # tool instead. It is slower (due to slow ArcGIS cursors) but
865        # does not crash.
866
867        from GeoEco.DataManagement.Fields import Field
868        Field.CalculateFields(arcGISConn, os.path.join(outputGDB, u'DensityRasters'), [u'StartTime', u'EndTime'], [u"(datetime.datetime(int(row.Name[1:5]), 1, 1, int(row.Name[8:10]), int(row.Name[10:12])) + datetime.timedelta(days=int(row.Name[5:8])-1)).strftime('%Y-%m-%d %H:%M:%S')", u"(datetime.datetime(int(row.Name[1:5]), 1, 1, int(row.Name[8:10]), int(row.Name[10:12])) + datetime.timedelta(days=int(row.Name[5:8])-1) + datetime.timedelta(seconds=" + unicode(simulationTimeStep*3600.0*summarizationPeriod) + u")).strftime('%Y-%m-%d %H:%M:%S')"], modulesToImport=[u'datetime'])
869
870        # Return successfully.
871
872        return outputDirectory
873
874###############################################################################
875# Metadata: module
876###############################################################################
877
878from GeoEco.ArcGIS import ArcGISDependency, ArcGISExtensionDependency
879from GeoEco.Dependencies import PythonAggregatedModuleDependency
880from GeoEco.Matlab import MatlabDependency
881from GeoEco.Metadata import *
882from GeoEco.Types import *
883
884AddModuleMetadata(shortDescription=_(u'Implements Eric Treml\'s coral reef connectivity analysis.'))
885
886###############################################################################
887# Metadata: CoralReefConnectivity class
888###############################################################################
889
890AddClassMetadata(CoralReefConnectivity,
891    shortDescription=_(u'Implements Eric Treml\'s coral reef connectivity analysis.'),
892    isExposedAsCOMServer=True,
893    comIID=u'{1B3A8B7F-8921-4A69-AF39-D209FAA0D588}',
894    comCLSID=u'{BAD2EBA8-1077-4A5B-B737-D4906466B60A}')
895
896# Public method: CoralReefConnectivity.CreateSimulationFromArcGISRasters
897
898AddMethodMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters,
899    shortDescription=_(u'Creates a coral reef connectivity simulation and initializes it with reef data in ArcGIS rasters.'),
900    isExposedToPythonCallers=True,
901    isExposedByCOM=True,
902    isExposedAsArcGISTool=True,
903    arcGISDisplayName=_(u'Create Coral Reef Connectivity Simulation From ArcGIS Rasters'),
904    arcGISToolCategory=_(u'Connectivity Analysis\\Coral Reef Larval Connectivity'),
905    dependencies=[ArcGISDependency(9, 2), ArcGISExtensionDependency(u'spatial'), PythonAggregatedModuleDependency('numpy')])
906
907AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'cls',
908    typeMetadata=ClassOrClassInstanceTypeMetadata(cls=CoralReefConnectivity),
909    description=_(u'%s class or an instance of it.') % CoralReefConnectivity.__name__)
910
911AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'simulationDirectory',
912    typeMetadata=DirectoryTypeMetadata(deleteIfParameterIsTrue=u'overwriteExisting', createParentDirectories=True),
913    description=_(
914u"""Output directory to create to contain the simulation's data.
915
916After creating the simulation directory, you must load ocean currents
917data into it using other tools before you can run the simulation. Use
918the MGET tools designed for this purpose. Unless you know what you are
919doing, do not modify the contents of the simulation directory
920yourself."""),
921    direction=u'Output',
922    arcGISDisplayName=_(u'Simulation directory to create'))
923
924AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'reefIDsRaster',
925    typeMetadata=ArcGISRasterLayerTypeMetadata(mustExist=True, allowedPixelTypes=[u'S8', u'U8', u'S16', u'U16', u'S32', u'U32', u'S64', u'U64']),
926    description=_(
927u"""Integer raster specifying the locations and IDs of reefs.
928
929The raster value specifies the ID of the reef contained by that cell.
930NoData indicates that the cell does not contain a reef. The value 0
931may not be used as a reef ID.
932
933The raster also defines the coordinate system, extent, and cell size
934for the analysis. The raster must use a projected coordinate system,
935with meters as the linear unit. The reef cover and water mask rasters
936must have the same coordinate system, extent, and cell size. When
937ocean currents data are loaded into the simulation, they will be
938automatically projected and clipped as needed into this coordinate
939system, extent, and cell size."""),
940    arcGISDisplayName=_(u'Reef IDs raster'))
941
942AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'reefCoverRaster',
943    typeMetadata=ArcGISRasterLayerTypeMetadata(mustExist=True),
944    description=_(
945u"""Floating-point raster specifying the proportion of each cell's
946area that is occupied by reef.
947
948This raster determines the density of coral larvae released into the
949cell at the start of the simulation. The density for a cell is
950determined by multiplying this raster's value by the larvae density
951parameter you specify when you run the simulation.
952
953This raster must have the same coordinate system, extent, and cell
954size as the reef IDs raster.
955
956The raster values must be greater than or equal to 0 and less than or
957equal to 1. The value 1 indicates that the entire cell is occupied by
958reef, while 0.5 indicates that only half of it is. If the cell size
959was 25 km by 25 km, this would mean the cell contained either 625 or
960312.5 square km of reef, respectively.
961
962If the value is 0 or NoData, it is assumed that the cell does not
963contain any reef, even if the reef IDs raster indicates that a reef is
964there. If the value is greater than 0 but the corresponding reef IDs
965raster contains NoData, it is assumed that the cell does not contain
966any reef."""),
967    arcGISDisplayName=_(u'Reef cover raster'))
968
969AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'waterMaskRaster',
970    typeMetadata=ArcGISRasterLayerTypeMetadata(mustExist=True),
971    description=_(
972u"""Raster specifying the locations of land and water. 0 or NoData
973indicates land; any other value indicates water.
974
975This raster must have the same coordinate system, dimensions, and cell
976size as the reef IDs raster."""),
977    arcGISDisplayName=_(u'Water mask raster'))
978
979AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'crosses180',
980    typeMetadata=BooleanTypeMetadata(),
981    description=_(
982u"""Set this option to True if your study area crosses the 180th
983meridian (i.e. 180 W / 180 E). This will happen if you are studying
984reefs in the tropical Pacific, for example."""),
985    arcGISDisplayName=_(u'Study area crosses the 180th meridian'))
986
987AddArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'overwriteExisting',
988    typeMetadata=BooleanTypeMetadata(),
989    description=_(
990u"""If True, the simulation directory will be deleted and recreated,
991if it exists. If False, a ValueError will be raised if the simulation
992directory exists."""),
993    initializeToArcGISGeoprocessorVariable=u'OverwriteOutput')
994
995# Public method: CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation
996
997AddMethodMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation,
998    shortDescription=_(u'Downloads Aviso geostrophic currents into a coral reef connectivity simulation.'),
999    longDescription=_AvisoGriddedProduct_LongDescription % {u'name': 'tool'},
1000    isExposedToPythonCallers=True,
1001    isExposedByCOM=True,
1002    isExposedAsArcGISTool=True,
1003    arcGISDisplayName=_(u'Load Aviso Geostrophic Currents Into Coral Reef Connectivity Simulation'),
1004    arcGISToolCategory=_(u'Connectivity Analysis\\Coral Reef Larval Connectivity'),
1005    dependencies=[ArcGISDependency(9, 2), ArcGISExtensionDependency(u'spatial'), PythonAggregatedModuleDependency('numpy')])
1006
1007CopyArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'cls', CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'cls')
1008
1009AddArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'simulationDirectory',
1010    typeMetadata=DirectoryTypeMetadata(mustExist=True),
1011    description=_(
1012u"""Existing coral reef connectivity simulation directory that should
1013recieve the geostrophic currents data.
1014
1015The directory must have been created using the Create Coral Reef
1016Connectivity Simulation tool.
1017
1018If you have already loaded geostrophic currents data into the
1019simulation, you can use this tool to add data for an additional range
1020of dates so long as the data is from the same Aviso product. If you
1021try to add a different product to existing data, this tool will report
1022an error.
1023
1024If you try to load data for a range of dates that have been already
1025loaded into the simulation, the downloads for those dates will be
1026skipped. There is no way to delete or overwrite the data that has
1027already been loaded into the simulation. If you want to overwrite your
1028existing data (e.g. you want to switch to using a different Aviso
1029product), you must create a new simulation."""),
1030    arcGISDisplayName=_(u'Simulation directory to recieve geostrophic currents data'))
1031
1032CopyArgumentMetadata(AvisoGriddedGeostrophicCurrents.CreateArcGISRasters, u'username', CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'username')
1033CopyArgumentMetadata(AvisoGriddedGeostrophicCurrents.CreateArcGISRasters, u'password', CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'password')
1034
1035AddArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'product',
1036    typeMetadata=UnicodeStringTypeMetadata(allowedValues=filter(lambda s: s.find(u'MADT') >= 0, AvisoGriddedGeostrophicCurrents.Products)),
1037    description=_(
1038u"""Aviso geostrophic currents product to load into the simulation.
1039Please see the Aviso web site (http://www.aviso.oceanobs.com) for
1040descriptions of the products.
1041
1042At the time of this writing, this tool supported all of the
1043geostrophic currents products that Aviso made available for download
1044over OPeNDAP. If you find that Aviso publishes a new product that is
1045not available using this tool, please contact the author of this tool
1046to have support for it added.
1047
1048The product name that you must pass for this parameter is case
1049sensitive. If you invoke this tool programmatically, be sure to
1050specify the product name using the proper case."""),
1051    arcGISDisplayName=_(u'Aviso geostrophic currents product to download'))
1052
1053AddArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'startDate',
1054    typeMetadata=DateTimeTypeMetadata(),
1055    description=_(
1056u"""Start date for the data to load into the simulation. Data that
1057occurs on or after the start date and on or before the end date will
1058be loaded into the simulation."""),
1059    arcGISDisplayName=_(u'Start date'))
1060
1061AddArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'endDate',
1062    typeMetadata=DateTimeTypeMetadata(),
1063    description=_(
1064u"""End date for the data to load into the simulation. Data that
1065occurs on or after the start date and on or before the end date will
1066be loaded into the simulation."""),
1067    arcGISDisplayName=_(u'End date'))
1068
1069CopyArgumentMetadata(AvisoGriddedGeostrophicCurrents.CreateArcGISRasters, u'timeout', CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'timeout')
1070CopyArgumentMetadata(AvisoGriddedGeostrophicCurrents.CreateArcGISRasters, u'maxRetryTime', CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'maxRetryTime')
1071CopyArgumentMetadata(AvisoGriddedGeostrophicCurrents.CreateArcGISRasters, u'cacheDirectory', CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'cacheDirectory')
1072
1073AddResultMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'updatedSimulationDirectory',
1074    typeMetadata=DirectoryTypeMetadata(),
1075    description=_(u'Updated simulation directory.'),
1076    arcGISDisplayName=_(u'Updated simulation directory'))
1077
1078# Public method: CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation
1079
1080AddMethodMetadata(CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation,
1081    shortDescription=_(u'Downloads HYCOM GOMl0.04 ocean currents into a coral reef connectivity simulation.'),
1082    longDescription=_HYCOMGOMl004_LongDescription % {u'name': 'tool'},
1083    isExposedToPythonCallers=True,
1084    isExposedByCOM=True,
1085    isExposedAsArcGISTool=True,
1086    arcGISDisplayName=_(u'Load HYCOM GOMl0.04 Currents Into Coral Reef Connectivity Simulation'),
1087    arcGISToolCategory=_(u'Connectivity Analysis\\Coral Reef Larval Connectivity'),
1088    dependencies=[ArcGISDependency(9, 2), ArcGISExtensionDependency(u'spatial'), PythonAggregatedModuleDependency('numpy')])
1089
1090CopyArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'cls', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'cls')
1091
1092AddArgumentMetadata(CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'simulationDirectory',
1093    typeMetadata=DirectoryTypeMetadata(mustExist=True),
1094    description=_(
1095u"""Existing coral reef connectivity simulation directory that should
1096recieve the currents.
1097
1098The directory must have been created using the Create Coral Reef
1099Connectivity Simulation tool.
1100
1101If you have already loaded HYCOM GOMl0.04 currents into the
1102simulation, you can use this tool to add currents for an additional
1103range of dates so long as they are for the same depth as the original
1104currents. If you try to add a different currents product to the
1105simulation (e.g. geostrophic currents from Aviso), or for a different
1106depth, this tool will report an error.
1107
1108If you try to load data for a range of dates that have been already
1109loaded into the simulation, the downloads for those dates will be
1110skipped. There is no way to delete or overwrite the data that has
1111already been loaded into the simulation. If you want to overwrite your
1112existing data (e.g. you want to switch to using a different depth),
1113you must create a new simulation."""),
1114    arcGISDisplayName=_(u'Simulation directory to recieve HYCOM currents'))
1115
1116CopyArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'startDate', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'startDate')
1117CopyArgumentMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'endDate', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'endDate')
1118
1119AddArgumentMetadata(CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'depth',
1120    typeMetadata=FloatTypeMetadata(allowedValues=[0., 5., 10., 15., 20., 25., 30., 40., 50., 60., 70., 80., 90., 100., 125., 150., 200., 250., 300., 400., 500., 600., 700., 800., 900., 1000., 1100., 1200., 1300., 1400., 1500., 1750., 2000., 2500., 3000., 3500., 4000., 4500., 5000., 5500.]),
1121    description=_(
1122u"""HYCOM depth layer to download currents from.
1123
1124This tool was designed primarily to study larvae that float at or near
1125the surface. The default depth is 0. If you are studying larvae that
1126stay submerged, you can choose a deeper depth, but be aware of two
1127important points:
1128
1129* This tool assumes the larvae remain at that depth for the entire
1130  simulation. It does not implement vertical migration or other
1131  behaviors that might be appropriate for your species. A newer
1132  version of this tool, to be released in 2011, will support vertical
1133  migration.
1134
1135* Be aware that the spatial extent of available currents data will
1136  shrink as depth increases. For example, if you choose a deep depth
1137  such as 250 m, there will be no data available for regions close to
1138  shore because the ocean is typically shallower than 250 m in those
1139  regions. Data for very deep depths will only be available near the
1140  center of the Gulf of Mexico.
1141"""),
1142    arcGISDisplayName=_(u'Minimum depth'))
1143
1144CopyArgumentMetadata(HYCOMGOMl0044D.CreateArcGISRasters, u'timeout', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'timeout')
1145CopyArgumentMetadata(HYCOMGOMl0044D.CreateArcGISRasters, u'maxRetryTime', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'maxRetryTime')
1146CopyArgumentMetadata(HYCOMGOMl0044D.CreateArcGISRasters, u'cacheDirectory', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'cacheDirectory')
1147
1148CopyResultMetadata(CoralReefConnectivity.LoadAvisoGeostrophicCurrentsIntoSimulation, u'updatedSimulationDirectory', CoralReefConnectivity.LoadHYCOMGOMl0044DCurrentsIntoSimulation, u'updatedSimulationDirectory')
1149
1150# Public method: CoralReefConnectivity.RunSimulation
1151
1152AddMethodMetadata(CoralReefConnectivity.RunSimulation,
1153    shortDescription=_(u'Executes a coral coral reef connectivity simulation.'),
1154    isExposedToPythonCallers=True,
1155    isExposedByCOM=True,
1156    isExposedAsArcGISTool=True,
1157    arcGISDisplayName=_(u'Run Coral Reef Connectivity Simulation'),
1158    arcGISToolCategory=_(u'Connectivity Analysis\\Coral Reef Larval Connectivity'),
1159    dependencies=[ArcGISDependency(9, 2), ArcGISExtensionDependency(u'spatial'), PythonAggregatedModuleDependency('numpy')])
1160
1161CopyArgumentMetadata(CoralReefConnectivity.CreateSimulationFromArcGISRasters, u'cls', CoralReefConnectivity.RunSimulation, u'cls')
1162
1163AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'simulationDirectory',
1164    typeMetadata=DirectoryTypeMetadata(mustExist=True),
1165    description=_(
1166u"""Existing coral reef connectivity simulation directory that has
1167been loaded with ocean currents data.
1168
1169The directory must have been created using the Create Coral Reef
1170Connectivity Simulation tool and then loaded with ocean currents data
1171using one of the tools provided for this purpose.
1172
1173The simulation may take hours to complete, due to the complex
1174mathematics involved. The run time is goverened by the size of your
1175study area, the number of reefs include in the simulation, and the
1176duration and time step of the simulation."""),
1177    arcGISDisplayName=_(u'Simulation directory'))
1178
1179AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'outputDirectory',
1180    typeMetadata=DirectoryTypeMetadata(mustExist=True),
1181    description=_(
1182u"""Output directory to receive the results of the simulation.
1183
1184The directory must already exist. Within it, the simulator creates:
1185
1186* DensityRasters - this subdirectory will contain a time series of
1187  rasters that represent snapshots of the larvae density throughout
1188  the study area, in particles per square km. A raster will be created
1189  each time the simulation is summarized; the summarization frequency
1190  is controlled by the Simulation Summarization Period parameter. The
1191  rasters' names will be of the form dYYYYDDDHHMM, where YYYY is the
1192  year, DDD is the day of the year (e.g. February 1 is day 032), HH is
1193  the hour, and MM is the minute.
1194
1195* DensityRasters raster catalog in ConnectivityGeodatabase.mdb - this
1196  catalog references the rasters created in the DensityRasters
1197  subdirectory. The StartDate and EndDate columns of the catalog
1198  specify the applicable timeframe for each raster, and may be used to
1199  animate the raster time series using the Animation toolbar available
1200  in ArcMap and other ArcGIS applications.
1201
1202* Edges feature class in ConnectivityGeodatabase.mdb - this line
1203  feature class shows which reefs are connected by larval dispersal.
1204  Each line represents a direction link between two reefs. The
1205  FromReefID and ToReefID fields specify the IDs of the source reef
1206  and sink reef, respectively. The MaxDispersal field shows the
1207  maximum quantity of larvae from the source reef found over the sink
1208  reef during the series of simulation summaries. The quantity is
1209  expressed as the fraction of larvae initially released from the
1210  source reef that are over the sink reef, and ranges from 0.0 to 1.0.
1211  For example, the value 0.01 means that 1% of the larvae initially
1212  released by the source reef were found over the sink reef, when
1213  dispersal from the source to the sink peaked.
1214"""),
1215    arcGISDisplayName=_(u'Output directory'))
1216
1217AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'startDate',
1218    typeMetadata=DateTimeTypeMetadata(),
1219    description=_(
1220u"""Start date of the simulation.
1221
1222The larvae are released from the reefs on this date. Selecting an
1223appropriate value requires knowledge of the reproductive biology of
1224the species you are studying."""),
1225    arcGISDisplayName=_(u'Start date'))
1226
1227AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'duration',
1228    typeMetadata=FloatTypeMetadata(minValue=0.0),
1229    description=_(
1230u"""Duration of the simulation, in days.
1231
1232Larger values require more computer memory and require more processing
1233time. To avoid extreme scenarios, this parameter can be no larger that
1234365 days."""),
1235    arcGISDisplayName=_(u'Duration'))
1236
1237AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'simulationTimeStep',
1238    typeMetadata=FloatTypeMetadata(minValue=0.0),
1239    description=_(
1240u"""Time step of the simulation, in hours.
1241
1242The time step defines the simulated time period at which larvae
1243density is recalculated using a Eulerian advection/diffusion model.
1244Smaller time steps increase the stability of the model and accuracy of
1245the results but also the run time and computer memory requirements of
1246the simulation. The original study from which this tool was developed
1247(Treml et al. 2008) used a time step of 2.4 hours.
1248
1249To check the model stability, this tool reports the effective Courant
1250number for the model, calculated from the simulation time step, grid
1251cell size, and average current velocity. The Courant number reflects
1252the portion of a cell that the larvae will traverse by advection in
1253one time step. If the Courant number is greater than 1, the model is
1254likely to be unstable and inaccurate. If it is significantly less than
12551 (e.g. 0.1), it is likely to be stable and accurate."""),
1256    arcGISDisplayName=_(u'Temp step'),
1257    arcGISCategory=_(u'Simulation options'))
1258
1259AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'summarizationPeriod',
1260    typeMetadata=IntegerTypeMetadata(minValue=1),
1261    description=_(
1262u"""Period, expressed as a number of simulation time steps, at which
1263the simulation should be summarized.
1264
1265The period specifies how frequently summarization should occur. The
1266first summarization will occur at the start of the simulation, and
1267subsequent summarizations will occur every time the period elapses.
1268For example, if the time step is 1 hour and the summarization period
1269is 24, summarizations will occur every 24 hours.
1270
1271The summarization procedure governs the production of simulation
1272outputs. Please see the documentation of the Output Directory
1273parameter for more information."""),
1274    arcGISDisplayName=_(u'Sumulation summarization period'),
1275    arcGISCategory=_(u'Simulation options'))
1276
1277AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'initialLarvaeDensity',
1278    typeMetadata=FloatTypeMetadata(minValue=1.0),
1279    description=_(
1280u"""Density of larvae, in particles per square km, released at the
1281start of the simulation.
1282
1283The original study from which this tool was developed (Treml et al.
12842008) used 10,000 particles per square km. This is consistent with
1285previous studies (Cowen et al. 2000; James et al. 2002; Largier 2003)
1286and was selected based on exploratory modeling efforts using densities
1287ranging from 1,000 to 100,000 particles per square km. For particular
1288taxa, these initial densities may need to be larger or smaller, based
1289on the specific fecundity and density characteristics of that species
1290of interest (Richmond 1987; Largier 2003)."""),
1291    arcGISDisplayName=_(u'Initial larvae density'),
1292    arcGISCategory=_(u'Simulation options'))
1293
1294AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'densityRasterCutoff',
1295    typeMetadata=FloatTypeMetadata(minValue=0.0, maxValue=100.0, canBeNone=True),
1296    description=_(
1297u"""Minimum larvae density, expressed as a percentage of the initial
1298larvae density, that a larvae density raster cell must be to not be
1299masked.
1300
1301This parameter is a percent of the initial larvae density. For
1302example, if the initial larvae density is 10000 particles per km and
1303this parameter is 1, a density of 100 particles per km will be used as
1304the cutoff. Larvae density cells that are less than this value will
1305set to NoData.
1306
1307This parameter does not affect the computations performed during the
1308simulation. It only affects the production of density rasters after
1309the simulation is complete, and is provided as a convenience for
1310visualization.
1311
1312If you set this parameter to 0, the density rasters will be produced
1313with no masking. The results may seem surprising. For most
1314simulations, larvae will have spread throughout the entire study area,
1315albeit in very small quantities for most cells. This is a surprising
1316effect of the diffusion component of the Eularian hydrodynamic
1317calculations. Diffiusion occurs equally in all directions at the rate
1318specified by the Diffusivity parameter. Given enough time, an
1319infinitesimal fraction of larvae from any given reef can theoretically
1320spread throughout the entire ocean simply by diffusion. To avoid a
1321confusing visualization, use this parameter to mask the extremely
1322density values that result from diffiusion."""),
1323    arcGISDisplayName=_(u'Cutoff percentage for larvae density rasters'),
1324    arcGISCategory=_(u'Simulation options'))
1325
1326##AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'deathRate',
1327##    typeMetadata=FloatTypeMetadata(minValue=0.0, maxValue=1.0, canBeNone=True),
1328##    description=_(
1329##u"""Death rate of the larvae, as proportion of the population that
1330##dies per day.
1331##
1332##Reliable mortality estimates for invertebrate larvae are rare due to
1333##the difficulty in sampling and identifying the same larval cohort in
1334##the plankton through time (Rumrill 1990; Morgan 1995). In a review of
1335##larval mortality rates, Rumrill (1990) showed that this parameter
1336##varied greatly between invertebrate species, from 0.016 to .357 per
1337##day, with an average of .223 per day. A recent study on several
1338##Pacific corals measuring larval survival, recruitment rates, and gene
1339##flow, reports the mortality of a broadcast-spawning coral to be around
1340##4 - 6% per day (Nishikawa et al. 2003; Nishikawa and Sakai 2005).  For
1341##this reason, the study from which this tool was developed (Treml et
1342##al. 2008) used a constant mortality of 6% per day. Although this rate
1343##is lower than the 18% per day reported for reef fish (Cowen et al.
1344##2000; James et al. 2002), it may be more representative of
1345##invertebrate taxa (Rumrill 1990, Ellien et al. 2004)."""),
1346##    arcGISDisplayName=_(u'Larvae death rate'),
1347##    arcGISCategory=_(u'Simulation options'))
1348
1349AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'diffusivity',
1350    typeMetadata=FloatTypeMetadata(minValue=0.0),
1351    description=_(
1352u"""The diffusion coefficent, in meters squared per second, to use in
1353the simulation.
1354
1355It is recommended that you consult an oceanographer to determine this
1356value. The original study from which this tool was developed (Treml et
1357al. 2008) used the value 25."""),
1358    arcGISDisplayName=_(u'Diffusivity'),
1359    arcGISCategory=_(u'Simulation options'))
1360
1361AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'includeReefIDs',
1362    typeMetadata=ListTypeMetadata(elementType=IntegerTypeMetadata(), canBeNone=True, minLength=1),
1363    description=_(
1364u"""List of IDs of the reefs to include in the simulation. If the list
1365is empty, all of the reefs will be included, unless some reefs are
1366listed in the "exclude reefs" parameter."""),
1367    arcGISDisplayName=_(u'IDs of reefs to include in the simulation'),
1368    arcGISCategory=_(u'Simulation options'))
1369
1370AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'excludeReefIDs',
1371    typeMetadata=ListTypeMetadata(elementType=IntegerTypeMetadata(), canBeNone=True, minLength=1),
1372    description=_(
1373u"""List of IDs of the reefs to exclude from the simulation. This list
1374may only be provided if the "include reefs" list is omitted (an error
1375will be reported if both lists are provided). If this list is
1376provided, all of the reefs will be included in the simulation except
1377those in this list."""),
1378    arcGISDisplayName=_(u'IDs of reefs to exclude from the simulation'),
1379    arcGISCategory=_(u'Simulation options'))
1380
1381AddArgumentMetadata(CoralReefConnectivity.RunSimulation, u'overwriteExisting',
1382    typeMetadata=BooleanTypeMetadata(),
1383    description=_(
1384u"""If True, the outputs will be overwritten, if it exists. If False,
1385a ValueError will be raised if the outputs exist."""),
1386    initializeToArcGISGeoprocessorVariable=u'OverwriteOutput')
1387
1388AddResultMetadata(CoralReefConnectivity.RunSimulation, u'updatedOutputDirectory',
1389    typeMetadata=DirectoryTypeMetadata(),
1390    description=_(u'Updated output directory.'),
1391    arcGISDisplayName=_(u'Updated output directory'))
1392
1393###############################################################################
1394# Names exported by this module
1395###############################################################################
1396
1397__all__ = ['CoralReefConnectivity']
Note: See TracBrowser for help on using the browser.