PopulateModel3D.java

Go to the documentation of this file.

//- ****************************************************************************
//- 
//- Copyright 2009 Sandia Corporation. Under the terms of Contract
//- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
//- retains certain rights in this software.
//- 
//- BSD Open Source License.
//- All rights reserved.
//- 
//- Redistribution and use in source and binary forms, with or without
//- modification, are permitted provided that the following conditions are met:
//- 
//-    * Redistributions of source code must retain the above copyright notice,
//-      this list of conditions and the following disclaimer.
//-    * Redistributions in binary form must reproduce the above copyright
//-      notice, this list of conditions and the following disclaimer in the
//-      documentation and/or other materials provided with the distribution.
//-    * Neither the name of Sandia National Laboratories nor the names of its
//-      contributors may be used to endorse or promote products derived from
//-      this software without specific prior written permission.
//- 
//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//- POSSIBILITY OF SUCH DAMAGE.
//-
//- ****************************************************************************
 
package gov.sandia.geotess.examples;
 
import gov.sandia.geotess.GeoTessMetaData;
import gov.sandia.geotess.GeoTessModel;
import gov.sandia.geotess.GeoTessModelUtils;
import gov.sandia.gmp.util.globals.DataType;
import gov.sandia.gmp.util.numerical.vector.VectorGeo;
 
import java.io.File;
import java.util.Date;
 
/**
 * An example of how to generate a 3D GeoTessModel and populate it with data.
 * The data are stored on a GeoTessGrid comprised of 3 multi-level 
 * tessellations, one for the core, one for the mantle and one for the crust.
 * The GeoTessGrid is not generated by this example.  It was previously 
 * computed using the GeoTessBuilder software and delivered as part of the 
 * GeoTess package in the GeoTessModels directory.
 * <p>The data used to populate the model come from the ak135 model,
 * which is hardcoded into the source code for the example.
 * 
 * @author sballar
 * 
 */
public class PopulateModel3D
{
/**
 * An example of how to generate a 3D GeoTessModel and populate it with data.
 * The data are stored on a GeoTessGrid comprised of 3 multi-level 
 * tessellations, one for the core, one for the mantle and one for the crust.
 * The GeoTessGrid is not generated by this example.  It was previously 
 * computed using the GeoTessBuilder software and delivered as part of the 
 * GeoTess package in the GeoTessModels directory.
 * <p>The data used to populate the model come from the ak135 model,
 * which is hardcoded into the source code for the example.
  * 
  * @param args this example requires a single command line argument
  * that specifies the full path of the grid file.  The file is called
  * small_model_grid.ascii and was delivered with the GeoTess package in 
  * the GeoTessModels directory.
  */
 public static void main(String[] args)
 {
  try
  {
   if (args.length == 0)
    throw new Exception(
      "\nMust specify a single command line argument specifying " +
      "the path to the file small_model_grid.geotess\n");
 
   System.out.println("Example that illustrates how to populate a 3D model.");
   System.out.println();
 
   // Create a MetaData object in which we can specify information
   // needed for model construction.
   GeoTessMetaData metaData = new GeoTessMetaData();
 
   // Specify a description of the model. This information is not
   // processed in any way by GeoTess. It is carried around for
   // information purposes.
   metaData.setDescription(String
     .format("Simple example of populating a 3D GeoTess model%n" +
       "comprised of 3 multi-level tessellations%n"
       + "author: Sandy Ballard%n"
       + "contact: sballar@sandia.gov"));
 
   // Specify a list of layer names delimited by semi-colons 
   metaData.setLayerNames("INNER_CORE; OUTER_CORE; LOWER_MANTLE; " +
     "TRANSITION_ZONE; UPPER_MANTLE; LOWER_CRUST; UPPER_CRUST");
 
   // Specify the relationship between grid tessellations and model layers.
   // the list has nLayers elements where each element specifies the 
   // index of the multilevel tessellation that supports the
   // corresponding layer.  In this example, the model has 
   // 7 layers and 3 multi-level tessellations.  The first 
   // tessellation (index 0) supports the core layers (layers 
   // 0 and 1).  The second tessellation (index 1) supports the 
   // three mantle layers (layer indices 2, 3 and 4).  The third
   // and final tessellation (index 2) supports the two crustal
   // layers (layer indices 5 and 6).
   metaData.setLayerTessIds(new int[] {0, 0, 1, 1, 1, 2, 2});
 
   // specify the names of the attributes and the units of the
   // attributes in two Strings delimited by semi-colons. 
   metaData.setAttributes("Vp; Vs; rho", "km/sec; km/sec; g/cc");
 
   // specify the DataType for the data. All attributes, in all
   // profiles, will have the same data type.  Note that this 
   // applies only to the data; radii are always stored as floats.
   metaData.setDataType(DataType.FLOAT);
 
   // specify the name of the software that is going to generate
   // the model.  This gets stored in the model for future reference.
   metaData.setModelSoftwareVersion("PopulateModel3D 1.0.0");
 
   // specify the date when the model was generated.  This gets 
   // stored in the model for future reference.
   metaData.setModelGenerationDate(new Date().toString());
 
   // specify the path to the file containing the grid to be used for
   // this test.  This information was passed in as a command line
   // argument.  Grids were included in the software delivery and
   // are available from the GeoTess website.  Grids can also be
   // constructed using the GeoTessBuilder software.  The grid 
   // required for this example is in a file called 
   // small_model_grid.ascii which is in the GeoTessModels 
   // directory delivered with GeoTess.
   File gridFile = new File(args[0]);
 
   // call a GeoTessModel constructor to build the model. This will
   // load the grid, and initialize all the data structures to null.
   // To be useful, we will have to populate the data structures.
   GeoTessModel model = new GeoTessModel(gridFile.getCanonicalPath(), metaData);
   
   // Now we will populate the model with Profiles.  At this point, the 
   // model has a 2D array of Profile objects with dimensions
   // nVertices x nLayers with all the elements of the array initialized
   // to null.  We will now populate the Profiles array.
   //
   // now loop over every vertex in the grid, connected and not-connected.
   for (int vtx = 0; vtx < model.getNVertices(); ++vtx)
   {
    // retrieve the unit vector corresponding to the i'th vertex 
    double[] vertex = model.getGrid().getVertex(vtx);
 
    // find the latitude and longitude of vertex, in degrees
    double lat = VectorGeo.getLatDegrees(vertex);
    double lon = VectorGeo.getLonDegrees(vertex);
 
    // loop over the 7 layers of the model (inner_core, outer_core, etc)
    for (int layer=0; layer<model.getNLayers(); ++layer)
    {
     // get the profile radii for the current lat, lon, and layer.
     // This is a 1D array of radius values where the first radius is
     // the radius at the bottom of the layer, the last radius is the
     // radius of the top of the layer and there can be as many in between
     // as desired.  This function will be different for every real-life
     // application.
     float[] radii = getRadii(lat, lon, layer);
 
     // get the profile data values for current lat, lon, layer.
     // This is a 2D array of values, nNodes by nAttributes.
     // In this example, nAttributes is 3 corresponding to 
     // attributes vp, vs and rho.  The number of nodes can
     // vary as a function of vertex and layer.
     float[][] rawData = getRawData(lat, lon, layer);
 
     model.setProfile(vtx, layer, radii, rawData);
    }
   }
 
 
 
   // At this point, we have a fully functional GeoTessModel object
   // that we can work with.
 
   // print a bunch of information about the model to the screen.
   System.out.println(model.toString());
 
   // print a profile for a single vertex, spanning all layers.
   System.out.println(GeoTessModelUtils.profileToString(model, 1));
 
   //model.writeModel(new File(gridFile.getParentFile(), "small_model.ascii"), "*");
   
   System.out.println("Done.");
  }
  catch (Exception ex)
  {
   ex.printStackTrace();
  }
 
 }
 
 /**
  * Returns a 1D profile of monotonically increasing radius values
  * that define the radial positions of nodes along a radial 
  * profile through a single layer in the model.
  * <p>
  * For this example, we will return the radius positions of the
  * nodes in the AK135 model, stretched a little bit so that 
  * the top of the model will coincide with the radius of the 
  * WGS84 ellipsoid instead of the ak135 value of 6371 km.
  * @param lat the latitude of the profile
  * @param lon the longitude of the profile
  * @param layer the index of the layer.
  * @return 1D array of radius values, in km.
  */
 protected static float[] getRadii(double lat, double lon, int layer)
 {
  // convert lat, lon in degrees to unit vector.
  double[] vertex = VectorGeo.getVectorDegrees(lat, lon);
 
  // find the radius of the WGS84 ellipsoid at the latitude of vertex.
  double earthRadius = VectorGeo.getEarthRadius(vertex);
 
  // find a stretching factor that will stretch the radius values so that they
  // span the range from zero at the center of the earth to the radius of the
  // WGS84 ellipsoid at the surface of the Earth.  This is not geophysically 
  // realistic, but sufficient for this simplistic example.
  float stretch = (float)(earthRadius / 6371.);
 
  float[] radii = new float[ak135[layer].length];
  for (int i=0; i<radii.length; ++i)
   radii[i] = ak135[layer][i][0] * stretch;
 
  return radii;
 }
 
 /**
  * Retrieve a 2D array of floats with nNodes x nAttributes elements.
  * The number of attributes is 3, for the vp, vs and rho.  nDataArrays
  * varies in the different layers.  For core and mantle layers, 
  * nNodes will equal the number of radii in the corresponding 
  * layers of the AK135 model. For the crustal layers, nNodes 
  * will be one, reflecting the fact that the attribute values are
  * constant in the crustal layers of the ak135 model.
  * <p>
  * In this example, the data returned are independent of latitude and
  * longitude since ak135 is a '1D' model, but this will not generally
  * be true for real 3D models.  
  * @param lat
  * @param lon
  * @param layer
  * @return
  */
 protected static float[][] getRawData(double lat, double lon, int layer)
 {
  // get a reference to the ak135 model values for the layer of interest.
  float[][] ak135Layer = ak135[layer];
  
  // nNodes is the number of radial positions defined in this layer in the 
  // ak135 model.
  int nNodes = ak135Layer.length;
 
  // data is the nNodes x nAttributes array that will be returned.
  float[][] data = null;
  
  if (nNodes == 2 && ak135Layer[0][1] == ak135Layer[1][1])
  {
   // this layer consists of only two nodes and the vp values are equal, 
   // which is true from the two crustal layers. Make a float[1][3] array 
   // containing the values from only the first node.
   data = new float[1][3];
   data[0][0] = ak135Layer[0][1];
   data[0][1] = ak135Layer[0][2];
   data[0][2] = ak135Layer[0][3];
  }
  else
  {
   // this layer has more than 2 nodes or maybe it has two 
   // nodes but the vp values are not equal.  Make an nNodes
   // by 3 array of data values.
   data = new float[nNodes][3];
   for (int i=0; i<nNodes; ++i)
   {
    data[i][0] = ak135Layer[i][1];
    data[i][1] = ak135Layer[i][2];
    data[i][2] = ak135Layer[i][3];
   }
  }
  return data;
 }
 
 /**
  * A 3D array of floats with nLayers x nNodes x nAttributes+1
  * elements.  The attributes are radius in km, 
  * vp in km/sec, vs in km/sec and density in g/cc.
  * There are 7 layers corresponding to the inner core, 
  * outer core, lower mantle, transition zone, upper mantle,
  * lower crust and upper crust.  The number of nodes in 
  * each layer is variable by layer.
  */
 static float[][][] ak135 = new float[][][] {{
  // inner core:
  // radius      vp          vs       density
  {   0.000F,  11.2622F,   3.6678F,  13.0122F},
  {  50.710F,  11.2618F,   3.6675F,  13.0117F},
  { 101.430F,  11.2606F,   3.6667F,  13.0100F},
  { 152.140F,  11.2586F,   3.6653F,  13.0074F},
  { 202.850F,  11.2557F,   3.6633F,  13.0036F},
  { 253.560F,  11.2521F,   3.6608F,  12.9988F},
  { 304.280F,  11.2477F,   3.6577F,  12.9929F},
  { 354.990F,  11.2424F,   3.6540F,  12.9859F},
  { 405.700F,  11.2364F,   3.6498F,  12.9779F},
  { 456.410F,  11.2295F,   3.6450F,  12.9688F},
  { 507.130F,  11.2219F,   3.6396F,  12.9586F},
  { 557.840F,  11.2134F,   3.6337F,  12.9474F},
  { 659.260F,  11.1941F,   3.6202F,  12.9217F},
  { 709.980F,  11.1832F,   3.6126F,  12.9072F},
  { 760.690F,  11.1715F,   3.6044F,  12.8917F},
  { 811.400F,  11.1590F,   3.5957F,  12.8751F},
  { 862.110F,  11.1457F,   3.5864F,  12.8574F},
  { 912.830F,  11.1316F,   3.5765F,  12.8387F},
  { 963.540F,  11.1166F,   3.5661F,  12.8188F},
  {1014.250F,  11.0983F,   3.5551F,  12.7980F},
  {1064.960F,  11.0850F,   3.5435F,  12.7760F},
  {1115.680F,  11.0718F,   3.5314F,  12.7530F},
  {1166.390F,  11.0585F,   3.5187F,  12.7289F},
  {1217.500F,  11.0427F,   3.5043F,  12.7037F}
 },
 {
  // outer core:
  // radius      vp          vs       density
  {1217.500F,  10.2890F,   0.0000F,  12.1391F},
  {1267.430F,  10.2854F,   0.0000F,  12.1133F},
  {1317.760F,  10.2745F,   0.0000F,  12.0867F},
  {1368.090F,  10.2565F,   0.0000F,  12.0593F},
  {1418.420F,  10.2329F,   0.0000F,  12.0311F},
  {1468.760F,  10.2049F,   0.0000F,  12.0001F},
  {1519.090F,  10.1739F,   0.0000F,  11.9722F},
  {1569.420F,  10.1415F,   0.0000F,  11.9414F},
  {1670.080F,  10.0768F,   0.0000F,  11.8772F},
  {1720.410F,  10.0439F,   0.0000F,  11.8437F},
  {1770.740F,  10.0103F,   0.0000F,  11.8092F},
  {1821.070F,   9.9761F,   0.0000F,  11.7737F},
  {1871.400F,   9.9410F,   0.0000F,  11.7373F},
  {1921.740F,   9.9051F,   0.0000F,  11.6998F},
  {1972.070F,   9.8682F,   0.0000F,  11.6612F},
  {2022.400F,   9.8304F,   0.0000F,  11.6216F},
  {2072.730F,   9.7914F,   0.0000F,  11.5809F},
  {2123.060F,   9.7513F,   0.0000F,  11.5391F},
  {2173.390F,   9.7100F,   0.0000F,  11.4962F},
  {2223.720F,   9.6673F,   0.0000F,  11.4521F},
  {2274.050F,   9.6232F,   0.0000F,  11.4069F},
  {2324.380F,   9.5777F,   0.0000F,  11.3604F},
  {2374.720F,   9.5306F,   0.0000F,  11.3127F},
  {2425.050F,   9.4814F,   0.0000F,  11.2639F},
  {2475.380F,   9.4297F,   0.0000F,  11.2137F},
  {2525.710F,   9.3760F,   0.0000F,  11.1623F},
  {2576.040F,   9.3205F,   0.0000F,  11.1095F},
  {2626.370F,   9.2634F,   0.0000F,  11.0555F},
  {2676.700F,   9.2042F,   0.0000F,  11.0001F},
  {2727.030F,   9.1426F,   0.0000F,  10.9434F},
  {2777.360F,   9.0792F,   0.0000F,  10.8852F},
  {2827.700F,   9.0138F,   0.0000F,  10.8257F},
  {2878.030F,   8.9461F,   0.0000F,  10.7647F},
  {2928.360F,   8.8761F,   0.0000F,  10.7023F},
  {2978.690F,   8.8036F,   0.0000F,  10.6385F},
  {3029.020F,   8.7283F,   0.0000F,  10.5731F},
  {3079.350F,   8.6496F,   0.0000F,  10.5062F},
  {3129.680F,   8.5692F,   0.0000F,  10.4378F},
  {3180.010F,   8.4861F,   0.0000F,  10.3679F},
  {3230.340F,   8.4001F,   0.0000F,  10.2964F},
  {3280.680F,   8.3122F,   0.0000F,  10.2233F},
  {3331.010F,   8.2213F,   0.0000F,  10.1485F},
  {3381.340F,   8.1283F,   0.0000F,  10.0722F},
  {3431.670F,   8.0382F,   0.0000F,   9.9942F},
  {3479.500F,   8.0000F,   0.0000F,   9.9145F}
 },
 {
  // lower mantle:
  // radius      vp          vs       density
  {3479.500F,  13.6602F,   7.2811F,   5.5515F},
  {3531.670F,  13.6566F,   7.2704F,   5.5284F},
  {3581.330F,  13.6530F,   7.2597F,   5.5051F},
  {3631.000F,  13.6494F,   7.2490F,   5.4817F},
  {3681.000F,  13.5900F,   7.2258F,   5.4582F},
  {3731.000F,  13.5312F,   7.2031F,   5.4345F},
  {3779.500F,  13.4741F,   7.1807F,   5.4108F},
  {3829.000F,  13.4156F,   7.1586F,   5.3869F},
  {3878.500F,  13.3585F,   7.1369F,   5.3628F},
  {3928.000F,  13.3018F,   7.1144F,   5.3386F},
  {3977.500F,  13.2465F,   7.0931F,   5.3142F},
  {4027.000F,  13.1894F,   7.0720F,   5.2898F},
  {4076.500F,  13.1336F,   7.0500F,   5.2651F},
  {4126.000F,  13.0783F,   7.0281F,   5.2403F},
  {4175.500F,  13.0222F,   7.0063F,   5.2154F},
  {4225.000F,  12.9668F,   6.9855F,   5.1904F},
  {4274.500F,  12.9096F,   6.9627F,   5.1652F},
  {4324.000F,  12.8526F,   6.9418F,   5.1398F},
  {4373.500F,  12.7956F,   6.9194F,   5.1143F},
  {4423.000F,  12.7382F,   6.8972F,   5.0887F},
  {4472.500F,  12.6804F,   6.8742F,   5.0629F},
  {4522.000F,  12.6221F,   6.8515F,   5.0370F},
  {4571.500F,  12.5631F,   6.8286F,   5.0109F},
  {4621.000F,  12.5031F,   6.8052F,   4.9847F},
  {4670.500F,  12.4426F,   6.7815F,   4.9584F},
  {4720.000F,  12.3819F,   6.7573F,   4.9319F},
  {4769.500F,  12.3185F,   6.7326F,   4.9052F},
  {4819.000F,  12.2550F,   6.7073F,   4.8785F},
  {4868.500F,  12.1912F,   6.6815F,   4.8515F},
  {4918.000F,  12.1245F,   6.6555F,   4.8245F},
  {4967.500F,  12.0577F,   6.6285F,   4.7973F},
  {5017.000F,  11.9895F,   6.6008F,   4.7699F},
  {5066.500F,  11.9200F,   6.5727F,   4.7424F},
  {5116.000F,  11.8491F,   6.5439F,   4.7148F},
  {5165.500F,  11.7766F,   6.5138F,   4.6870F},
  {5215.000F,  11.7026F,   6.4828F,   4.6591F},
  {5264.500F,  11.6269F,   6.4510F,   4.6310F},
  {5314.000F,  11.5495F,   6.4187F,   4.6028F},
  {5363.500F,  11.4705F,   6.3854F,   4.5744F},
  {5413.000F,  11.3896F,   6.3512F,   4.5459F},
  {5462.500F,  11.3068F,   6.3160F,   4.5173F},
  {5512.000F,  11.2221F,   6.2798F,   4.4885F},
  {5561.500F,  11.1353F,   6.2426F,   4.4596F},
  {5611.000F,  11.0558F,   6.2095F,   4.4305F},
  {5661.000F,  10.9229F,   6.0897F,   4.4010F},
  {5711.000F,  10.7900F,   5.9600F,   4.3714F}
 },
 {
  // transition zone:
  // radius      vp          vs       density
  {5711.000F,  10.2000F,   5.6100F,   4.0646F},
  {5761.000F,  10.0320F,   5.5040F,   4.0028F},
  {5811.000F,   9.8640F,   5.3980F,   3.9410F},
  {5861.000F,   9.6960F,   5.2920F,   3.8793F},
  {5911.000F,   9.5280F,   5.1860F,   3.8175F},
  {5961.000F,   9.3600F,   5.0800F,   3.7557F}
 },
 {
  // upper mantle:
  // radius      vp          vs       density
  {5961.000F,   9.0300F,   4.8700F,   3.5470F},
  {6011.000F,   8.8475F,   4.7830F,   3.5167F},
  {6061.000F,   8.6650F,   4.6960F,   3.4864F},
  {6111.000F,   8.4825F,   4.6090F,   3.4561F},
  {6161.000F,   8.3000F,   4.5230F,   3.4258F},
//  {6161.000F,   8.3000F,   4.5180F,   3.4258F}, // ignore small S discontinuity at 210 km depth
  {6206.000F,   8.1750F,   4.5090F,   3.3985F},
  {6251.000F,   8.0500F,   4.5000F,   3.3713F},
  {6293.500F,   8.0450F,   4.4900F,   3.3455F},
  {6336.000F,   8.0400F,   4.4800F,   3.3198F}
 },
 {
  // lower crust:
  // radius      vp          vs       density
  {6336.000F,   6.5000F,   3.8500F,   2.9200F},
  {6351.000F,   6.5000F,   3.8500F,   2.9200F}
 },
 {
  // upper crust:
  // radius      vp          vs       density
  {6351.000F,   5.8000F,   3.4600F,   2.7200F},
  {6371.000F,   5.8000F,   3.4600F,   2.7200F}}};
}