diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/BoundingHazards.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/BoundingHazards.java
new file mode 100644
index 0000000000000000000000000000000000000000..f89e2cf7898fb92411ac7209c331636b7d78ff50
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/BoundingHazards.java
@@ -0,0 +1,217 @@
+package gov.usgs.earthquake.nshmp.netcdf;
+
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.EnumMap;
+import java.util.List;
+import java.util.Map;
+
+import com.google.common.collect.Maps;
+
+import gov.usgs.earthquake.nshmp.data.XySequence;
+import gov.usgs.earthquake.nshmp.geo.Location;
+import gov.usgs.earthquake.nshmp.gmm.Imt;
+import gov.usgs.earthquake.nshmp.netcdf.NetcdfUtils.Key;
+
+import ucar.ma2.Array;
+import ucar.ma2.DataType;
+import ucar.ma2.InvalidRangeException;
+import ucar.nc2.Group;
+import ucar.nc2.Variable;
+import ucar.nc2.dataset.NetcdfDataset;
+
+/*
+ * Container for gridded hazard curves at four closest grid points to target
+ */
+class BoundingHazards {
+
+ final Netcdf netcdf;
+ final String baseGroup;
+ final NetcdfCoordinates coords;
+
+ private int idxLonLL;
+ private int idxLatLL;
+
+ final List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazards;
+
+ private BoundingHazards(Netcdf netcdf, Location site, String baseGroup) {
+ this.netcdf = netcdf;
+ this.baseGroup = baseGroup;
+ this.coords = netcdf.coordinates();
+
+ coords.checkCoords(site);
+
+ var longitudes = coords.locations().stream().mapToDouble(loc -> loc.longitude).toArray();
+ var latitudes = coords.locations().stream().mapToDouble(loc -> loc.latitude).toArray();
+
+ idxLonLL = NetcdfUtils.getIdxLTEQ(longitudes, site.longitude);
+ idxLatLL = NetcdfUtils.getIdxLTEQ(latitudes, site.latitude);
+
+ boundingHazards = extractHazardsAt(idxLonLL, idxLatLL);
+
+ var fracLon = NetcdfUtils.calcGridFrac(longitudes, idxLonLL, site.longitude);
+ var fracLat = NetcdfUtils.calcGridFrac(latitudes, idxLatLL, site.latitude);
+ boundingHazards.add(calcTargetHazards(fracLon, fracLat));
+ }
+
+ /**
+ * Returns the bounding hazards at four closet grid points to target.
+ *
+ * @param netcdf The {@code Netcdf}
+ * @param site The site to get bounding hazards
+ * @param baseGroup The Netcdf base group
+ */
+ static List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazards(
+ Netcdf netcdf,
+ Location site,
+ String baseGroup) {
+ return new BoundingHazards(netcdf, site, baseGroup).boundingHazards;
+ }
+
+ private Map<SiteClass, Map<Imt, XySequence>> calcTargetHazards(double fracLon, double fracLat) {
+ var westTarget = getTargetData(boundingHazards.get(0), boundingHazards.get(1), fracLat);
+ var eastTarget = getTargetData(boundingHazards.get(3), boundingHazards.get(2), fracLat);
+ return getTargetData(westTarget, eastTarget, fracLon);
+ }
+
+ /*
+ * Return hazard curves for four closest grid points as a List of:
+ * Map<SiteClass, Map<Imt, XySequence(iml,hazard)>>
+ *
+ * List order is clockwise from lower left corner: LL, UL, UR, LR, [T] with an
+ * empty slot for the interpolated target hazards
+ */
+ private List<Map<SiteClass, Map<Imt, XySequence>>> extractHazardsAt(int idxLonLL, int idxLatLL) {
+
+ var boundingHazardMaps = new ArrayList<Map<SiteClass, Map<Imt, XySequence>>>(5);
+
+ try (NetcdfDataset ncd = NetcdfDataset.openDataset(netcdf.path().toString())) {
+ Group targetGroup = ncd.findGroup(baseGroup);
+
+ // TODO: rename variable in netCDF
+ Variable vHazards = targetGroup.findVariable(Key.AEPS);
+
+ // set up origin and shape arrays for reading hazard curves at four
+ // bounding
+ // (lon,lat) grid points
+ int[] origin = new int[] { 0, 0, idxLatLL, idxLonLL, 0 };
+ int[] shape = new int[] {
+ coords.siteClasses().size(),
+ coords.imt().size(),
+ 2,
+ 2,
+ coords.nIml()
+ };
+
+ // read data into Array for bounding grid points
+ var aHazards = vHazards.read(origin, shape);
+
+ // Array aHazards now has shape [nVs,nImt,2,2,nIml]
+ // ...so origin will now be [0,0,0,0,0] for LL grid point
+ // ...and shape of requested array is [nVs,nImt,1,1,nIml]
+
+ // Array aHazards is size 2 on the lon and lat dimensions
+ shape[2] = 1;
+ shape[3] = 1;
+
+ // TODO: XySequence instead of double[]
+
+ // Put data into hazard map for LL point
+ origin[2] = 0; // lat index position
+ origin[3] = 0; // lon index position
+ boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
+ // Put data into hazard map for UL point
+ origin[2] = 1; // lat index position
+ origin[3] = 0; // lon index position
+ boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
+ // Put data into hazard map for UR point
+ origin[2] = 1; // lat index position
+ origin[3] = 1; // lon index position
+ boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
+ // Put data into hazard map for LR point
+ origin[2] = 0; // lat index position
+ origin[3] = 1; // lon index position
+
+ boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
+ } catch (IOException | InvalidRangeException e) {
+ throw new RuntimeException("Could not read Netcdf file [" + netcdf.path() + "]");
+ }
+
+ // interpolate hazards to target point
+ return boundingHazardMaps;
+
+ }
+
+ /*
+ * Read hazard curves from netCDF variable into map of hazards by SiteClass
+ * and Imt
+ *
+ * TODO: if target is on a grid point (or on a grid lat or lon), no need to
+ * read 4 bounding points ?
+ */
+ private Map<SiteClass, Map<Imt, XySequence>> mapHazardsFromArray(Array hazards) {
+ // hazards has had its lat and lon dimensions reduced (so it now is rank 3)
+ // hazards[SiteClass, Imt, Iml]
+
+ EnumMap<SiteClass, Map<Imt, XySequence>> vsImtHazardMap = Maps.newEnumMap(SiteClass.class);
+
+ for (int iSc = 0; iSc < coords.siteClasses().size(); iSc++) {
+ var sc = coords.siteClasses().get(iSc);
+
+ Map<Imt, XySequence> imtHazardMap = Maps.newEnumMap(Imt.class);
+ for (int iImt = 0; iImt < coords.imt().size(); iImt++) {
+ Imt imt = coords.imt().get(iImt);
+
+ // set origin and shape of double[] hazards to read
+ var origin = new int[] { iSc, iImt, 0 };
+ var shape = new int[] { 1, 1, coords.nIml() };
+
+ // try block needed for sectioning Array
+ try {
+ imtHazardMap.put(imt,
+ XySequence.create(
+ coords.iml().get(sc).get(imt),
+ (double[]) hazards.section(origin, shape).reduce()
+ .get1DJavaArray(DataType.DOUBLE)));
+ } catch (InvalidRangeException e) {
+ throw new RuntimeException(e.getMessage());
+ }
+ }
+
+ vsImtHazardMap.put(sc, imtHazardMap);
+ }
+
+ return Maps.immutableEnumMap(vsImtHazardMap);
+ }
+
+ /*
+ * Get data for target point
+ *
+ * @param d1 data at first point (p1)
+ *
+ * @param d2 data at second point (p2)
+ *
+ * @param frac fractional distance between p1 and p2 to target point
+ */
+ private static Map<SiteClass, Map<Imt, XySequence>> getTargetData(
+ Map<SiteClass, Map<Imt, XySequence>> d1,
+ Map<SiteClass, Map<Imt, XySequence>> d2,
+ double frac) {
+ // do we need better checking here? or is it safe to assume that every
+ // Map<SiteClass, Map<Imt,double[]>> passed in is consistent?
+ if (d1.size() != d2.size()) {
+ throw new IllegalArgumentException("Array size disagreement, cannot interpolate");
+ }
+
+ if (frac == 0.0) {
+ // target is the same as d1
+ return d1;
+ } else if (frac == 1.0) {
+ // target is the same as d2
+ return d2;
+ } else {
+ return NetcdfUtils.linearInterpolate(d1, d2, frac);
+ }
+ }
+
+}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/Netcdf.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/Netcdf.java
new file mode 100644
index 0000000000000000000000000000000000000000..048b8bcfd52e099363a5962f43a45bc6e3f31ab4
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/Netcdf.java
@@ -0,0 +1,60 @@
+package gov.usgs.earthquake.nshmp.netcdf;
+
+import java.io.IOException;
+import java.nio.file.Path;
+import java.util.List;
+import java.util.Map;
+
+import gov.usgs.earthquake.nshmp.data.XySequence;
+import gov.usgs.earthquake.nshmp.geo.Location;
+import gov.usgs.earthquake.nshmp.gmm.Imt;
+
+public interface Netcdf {
+
+ /**
+ * Returns a {@code NetCdf} reader for 2018 hazard curves.
+ *
+ * @param path The path to the 2018 Netcdf hazard file
+ * @throws IOException
+ */
+ public static Netcdf create2018Reader(Path path) {
+ return new Netcdf2018Reader(path);
+ }
+
+ /**
+ * Returns the Netcdf dimensions and coordinate variables.
+ */
+ NetcdfCoordinates coordinates();
+
+ /**
+ * Returns the path to the Netcdf file.
+ */
+ Path path();
+
+ /**
+ * Return the full set of hazard curves at the bounding grid points and the
+ * target Location. Intended for validation uses.
+ *
+ * @param site The location to get bouning hazards
+ */
+ List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazards(Location site);
+
+ /**
+ * Return a Map<SiteClass, Map<Imt, XySequence>> with hazard curves at the
+ * specified Location.
+ *
+ * @param site The site to get the hazard curves
+ */
+ Map<SiteClass, Map<Imt, XySequence>> hazard(Location site);
+
+ /**
+ * Return a single hazard curve for the specified Location, SiteClass, and
+ * Imt.
+ *
+ * @param site The site to get the hazard curves
+ * @param siteClass The site class
+ * @param imt The IMT
+ */
+ XySequence hazard(Location site, SiteClass siteClass, Imt imt);
+
+}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/Netcdf2018Reader.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/Netcdf2018Reader.java
new file mode 100644
index 0000000000000000000000000000000000000000..0e482fdb1bc193223d6c638896110e85b3fae7e2
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/Netcdf2018Reader.java
@@ -0,0 +1,87 @@
+package gov.usgs.earthquake.nshmp.netcdf;
+
+import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.util.List;
+import java.util.Map;
+
+import gov.usgs.earthquake.nshmp.data.XySequence;
+import gov.usgs.earthquake.nshmp.geo.Location;
+import gov.usgs.earthquake.nshmp.gmm.Imt;
+
+import ucar.nc2.dataset.NetcdfDataset;
+
+/**
+ *
+ * Hazard curves stored in netCDF
+ *
+ * This class is currently specific to NSHMP 2018a data release.
+ *
+ * TODO: generalize for any NSHMP model data release with subclasses per data
+ * release.
+ *
+ * TODO: Replace loops with streams
+ *
+ * TODO: Utilize XySequence.copyOf()
+ *
+ * TODO: Update netCDF file format, rename variables, implement enums,
+ * standardize attributes, remove Vs30 dimension from Imls
+ *
+ * TODO: attempt to calculate index, possibly as
+ * NshmpNetcdfCoordinates.calcIdxLTEQ()
+ *
+ * @author U.S. Geological Survey
+ */
+class Netcdf2018Reader implements Netcdf {
+
+ private static final String BASE_GROUP = "/CONUS/2018/v1.1";
+
+ private final Path path;
+ private final NetcdfCoordinates coords;
+
+ Netcdf2018Reader(Path path) {
+ this.path = path;
+
+ if (Files.notExists(path)) {
+ throw new IllegalArgumentException("Path to Netcdf file [" + path + "] does not exisit");
+ }
+
+ try (var ncd = NetcdfDataset.openDataset(path.toString())) {
+ // TODO:
+ // - Handle metadata from netCDF file get root group and
+ // attributes, and attributes of subgroups
+ // - Get netCDF dimensions
+ var targetGroup = ncd.findGroup(BASE_GROUP);
+ coords = new NetcdfCoordinates(targetGroup);
+ } catch (IOException e) {
+ throw new RuntimeException("Could not read Netcdf file [" + path + " ]");
+ }
+ }
+
+ @Override
+ public NetcdfCoordinates coordinates() {
+ return coords;
+ }
+
+ @Override
+ public Path path() {
+ return path;
+ }
+
+ @Override
+ public List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazards(Location site) {
+ return BoundingHazards.boundingHazards(this, site, BASE_GROUP);
+ }
+
+ @Override
+ public Map<SiteClass, Map<Imt, XySequence>> hazard(Location site) {
+ return boundingHazards(site).get(4);
+ }
+
+ @Override
+ public XySequence hazard(Location site, SiteClass siteClass, Imt imt) {
+ return hazard(site).get(siteClass).get(imt);
+ }
+
+}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NetcdfCoordinates.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NetcdfCoordinates.java
new file mode 100644
index 0000000000000000000000000000000000000000..43c6d3910479b17840b13dd4fe7ba4eee2fba764
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NetcdfCoordinates.java
@@ -0,0 +1,156 @@
+package gov.usgs.earthquake.nshmp.netcdf;
+
+import static com.google.common.base.Preconditions.checkArgument;
+
+import java.io.IOException;
+import java.util.Arrays;
+import java.util.EnumMap;
+import java.util.List;
+import java.util.Map;
+import java.util.stream.Collectors;
+
+import com.google.common.collect.Maps;
+
+import gov.usgs.earthquake.nshmp.geo.Location;
+import gov.usgs.earthquake.nshmp.geo.LocationList;
+import gov.usgs.earthquake.nshmp.gmm.Imt;
+import gov.usgs.earthquake.nshmp.netcdf.NetcdfUtils.Key;
+
+import ucar.ma2.DataType;
+import ucar.ma2.InvalidRangeException;
+import ucar.nc2.Group;
+import ucar.nc2.Variable;
+
+/*
+ * Container for dimensions and coordinate variables of nshmp netCDF file to
+ * facilitate indexing operations prior to data retrieval.
+ */
+public class NetcdfCoordinates {
+
+ private final List<SiteClass> siteClasses;
+ private final List<Imt> imts;
+ private final LocationList locations;
+ private final Map<SiteClass, Map<Imt, double[]>> imls;
+ private final int nIml;
+
+ NetcdfCoordinates(Group targetGroup) throws IOException {
+ // This bypasses the netCDF dimensions, but since we know what the
+ // variables and their dimensions should be, this is OK(???)
+
+ // get coordinate variables
+ var vVs30s = targetGroup.findVariable(Key.VS30);
+ var vImts = targetGroup.findVariable(Key.IMT);
+ var vlats = targetGroup.findVariable(Key.LAT);
+ var vlons = targetGroup.findVariable(Key.LON);
+ var vImls = targetGroup.findVariable(Key.IMLS);
+
+ // read coordinate variables into java arrays
+ var vs30s = (int[]) vVs30s.read().get1DJavaArray(DataType.INT);
+ var imtVals = (double[]) vImts.read().get1DJavaArray(DataType.DOUBLE);
+ var lats = (double[]) vlats.read().get1DJavaArray(DataType.DOUBLE);
+ var lons = (double[]) vlons.read().get1DJavaArray(DataType.DOUBLE);
+
+ var builder = LocationList.builder();
+ for (int i = 0; i < lons.length; i++) {
+ builder.add(lats[i], lons[i]);
+ }
+
+ locations = builder.build();
+
+ // vImls has dimensions (vs30, Imt, Iml)
+ // alternatively get nIml from Dimension Iml
+ nIml = vImls.getDimension(2).getLength(); // Imls.length;
+
+ // convert Vs30<int> to SiteClass Enum
+ // TODO: site class should be stored as enum in netCDF file rather than
+ // vs30 value, so this will break:
+ siteClasses = Arrays.stream(vs30s)
+ .mapToObj(vs30 -> SiteClass.ofValue(vs30))
+ .collect(Collectors.toUnmodifiableList());
+
+ // TODO: Imt should be stored in netCDF as enum, so this will break:
+ imts = Arrays.stream(imtVals)
+ .mapToObj(x -> (x < 0.009) ? Imt.PGA : Imt.fromPeriod(x))
+ .collect(Collectors.toUnmodifiableList());
+
+ imls = mapImls(vImls);
+ }
+
+ LocationList locations() {
+ return LocationList.copyOf(locations);
+ }
+
+ int nIml() {
+ return nIml;
+ }
+
+ List<SiteClass> siteClasses() {
+ return List.copyOf(siteClasses);
+ }
+
+ List<Imt> imt() {
+ return List.copyOf(imts);
+ }
+
+ Map<SiteClass, Map<Imt, double[]>> iml() {
+ return imls;
+ }
+
+ /*
+ * validate target coordinates
+ */
+ void checkCoords(Location site) {
+ var bounds = locations.bounds();
+ var minLon = bounds.min.longitude;
+ var maxLon = bounds.max.longitude;
+ var minLat = bounds.min.latitude;
+ var maxLat = bounds.max.latitude;
+ var lon = site.longitude;
+ var lat = site.latitude;
+
+ checkArgument(
+ (lon >= minLon && lon <= maxLon && lat >= minLat && lat <= maxLat),
+ String.format("Target coordinates out of range, %.3f <= lon < %.3f, %.3f <= lat <= %.3f",
+ minLon, maxLat, minLat, maxLat));
+ }
+
+ /*
+ * convert 3D Iml variable (dimensions vs30, Imt, Iml) to map of Imls by
+ * SiteClass and Imt
+ *
+ * TODO: use MultiMap or SetMultiMap (etc.) to store unique IML sets? Could
+ * then also initialize the underlying XySequence objects for reading in the
+ * hazard curves...
+ */
+ private Map<SiteClass, Map<Imt, double[]>> mapImls(Variable vImls) {
+
+ EnumMap<SiteClass, Map<Imt, double[]>> vsImtImlMap = Maps.newEnumMap(SiteClass.class);
+
+ for (int i = 0; i < siteClasses.size(); i++) {
+ var sc = siteClasses.get(i);
+
+ Map<Imt, double[]> imtImlMap = Maps.newEnumMap(Imt.class);
+ for (int j = 0; j < imts.size(); j++) {
+ var imt = imts.get(j);
+
+ // set origin and shape of double[] Imls to read
+ var origin = new int[] { i, j, 0 };
+ var shape = new int[] { 1, 1, nIml };
+
+ try {
+ imtImlMap.put(imt,
+ (double[]) vImls.read(origin, shape).reduce().get1DJavaArray(DataType.DOUBLE));
+ } catch (IOException | InvalidRangeException e) {
+ var msg = "Failed read attempt for vImls with origin: " +
+ Arrays.toString(origin) + ", shape: " + Arrays.toString(shape);
+ throw new RuntimeException(msg);
+ }
+ }
+
+ vsImtImlMap.put(sc, imtImlMap);
+ }
+
+ return Maps.immutableEnumMap(vsImtImlMap);
+ }
+
+}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NetcdfUtils.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NetcdfUtils.java
new file mode 100644
index 0000000000000000000000000000000000000000..28a712c6a78df6240272657a9b8687661eb6fcd6
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NetcdfUtils.java
@@ -0,0 +1,123 @@
+package gov.usgs.earthquake.nshmp.netcdf;
+
+import java.util.Arrays;
+import java.util.Map;
+
+import com.google.common.collect.Maps;
+import com.google.common.math.DoubleMath;
+
+import gov.usgs.earthquake.nshmp.data.XySequence;
+import gov.usgs.earthquake.nshmp.gmm.Imt;
+
+class NetcdfUtils {
+
+ /*
+ * find index of first element in a (sorted ascending) that is less than or
+ * equal to target value t. If target value is equal to the maximum value in a
+ * (the last element), an index of a.length - 1 is returned.
+ */
+ static int getIdxLTEQ(double[] a, double t) {
+ // assumes array is in sorted order (a[i] < a[i+1])
+ // make sure target is within the range of the array
+ var tol = 0.000001;
+ var n = a.length;
+
+ if (t < a[0] || t > a[n - 1]) {
+ // should never get here thanks to NshmpNetcdfCoordinates.checkCoords()
+ throw new IllegalArgumentException(
+ String.format("Target outside of valid range: %.4f <= t <= %4f", a[0], a[n - 1]));
+ }
+
+ // if (t == a[0]) {
+ if (DoubleMath.fuzzyEquals(a[0], t, tol)) {
+ return 0;
+ }
+ // if (t == a[n - 1]) {
+ if (DoubleMath.fuzzyEquals(a[n - 1], t, tol)) {
+ return n - 2; // return second to last index number
+ }
+
+ var idx = Arrays.binarySearch(a, t);
+ if (idx < 0) {
+ // "exact" match not found, so use index of first value less than target
+ // this is insertion_point - 1
+ // returned idx = -insertion_point - 1
+ idx = -1 * (idx + 1) - 1;
+ // array[idx] <= target
+ }
+
+ return idx;
+ }
+
+ /*
+ * Calculate fractional distance from a1 to t, between a1 and a2
+ */
+ static double calcFrac(double a1, double a2, double t) {
+ var tol = 0.00001;
+
+ if (Math.abs(t - a1) < tol) {
+ // target value == a1
+ return 0.0;
+ } else if (Math.abs(t - a2) < tol) {
+ // target value == a2
+ return 1.0;
+ } else {
+ // calculate fractional distance to t between a[i] and a[i+1]
+ return (t - a1) / (a2 - a1);
+ }
+ }
+
+ /*
+ * Calculate fractional distance from a[i] to t, between a[i] and a[i+1]
+ */
+ static double calcGridFrac(double[] a, int i, double t) {
+ return calcFrac(a[i], a[i + 1], t);
+ }
+
+ /*
+ * Linear interpolation of data values to a target point
+ *
+ * TODO: use XySequence.copyOf(), streams
+ */
+ static Map<SiteClass, Map<Imt, XySequence>> linearInterpolate(
+ Map<SiteClass, Map<Imt, XySequence>> v1,
+ Map<SiteClass, Map<Imt, XySequence>> v2,
+ double frac) {
+ // return v1 * (1 - frac) + v2 * frac
+ // DoubleData routines add and multiply in place, we don't want to modify
+ // the original data here
+ if (v1.size() != v2.size()) {
+ throw new IllegalArgumentException("Array size disagreement, cannot interpolate");
+ }
+
+ Map<SiteClass, Map<Imt, XySequence>> tMap = Maps.newEnumMap(SiteClass.class);
+
+ for (SiteClass sc : v1.keySet()) {
+ Map<Imt, XySequence> imtHazards = Maps.newEnumMap(Imt.class);
+
+ for (Imt imt : v1.get(sc).keySet()) {
+ var v1Haz = v1.get(sc).get(imt).yValues().toArray();
+ var v2Haz = v2.get(sc).get(imt).yValues().toArray();
+ var t = new double[v1Haz.length];
+
+ for (int i = 0; i < v1Haz.length; i++) {
+ t[i] = v1Haz[i] * (1 - frac) + v2Haz[i] * frac;
+ }
+ imtHazards.put(imt, XySequence.create(v1.get(sc).get(imt).xValues().toArray(), t));
+ }
+ tMap.put(sc, imtHazards);
+ }
+ return tMap;
+ }
+
+ static class Key {
+ static final String AEPS = "AEPs";
+ static final String IMLS = "Imls";
+ static final String IMT = "Imt";
+ static final String LAT = "lat";
+ static final String LON = "lon";
+ static final String POE = "Poe";
+ static final String SPACING = "spacing";
+ static final String VS30 = "vs30";
+ }
+}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NshmpNetcdfReader.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NshmpNetcdfReader.java
deleted file mode 100644
index 12e099c60f1075286a7407bec0f75680dcc1fe2e..0000000000000000000000000000000000000000
--- a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/NshmpNetcdfReader.java
+++ /dev/null
@@ -1,926 +0,0 @@
-package gov.usgs.earthquake.nshmp.netcdf;
-
-import static com.google.common.base.Preconditions.checkArgument;
-
-import java.io.IOException;
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.EnumMap;
-import java.util.List;
-import java.util.Map;
-import java.util.logging.Level;
-import java.util.logging.Logger;
-import java.util.stream.Collectors;
-
-import ucar.ma2.Array;
-import ucar.ma2.DataType;
-import ucar.ma2.InvalidRangeException;
-import ucar.nc2.Group;
-import ucar.nc2.Variable;
-import ucar.nc2.dataset.NetcdfDataset;
-
-import com.google.common.collect.Maps;
-import com.google.common.math.DoubleMath;
-
-import gov.usgs.earthquake.nshmp.data.XySequence;
-import gov.usgs.earthquake.nshmp.geo.Location;
-import gov.usgs.earthquake.nshmp.gmm.Imt;
-
-/**
- *
- * Interface to static hazard curves stored in netCDF
- *
- * This class is currently specific to NSHMP 2018a data release.
- *
- * TODO: generalize for any NSHMP model data release with subclasses per data
- * release.
- *
- * TODO: Replace loops with streams
- *
- * TODO: Utilize XySequence.copyOf()
- *
- * TODO: Update netCDF file format, rename variables, implement enums,
- * standardize attributes, remove Vs30 dimension from Imls
- *
- * @author U.S. Geological Survey
- */
-public class NshmpNetcdfReader {
-
- final String nshmpNetcdfFilename; // path to netCDF file
- final String nshmpBaseGroup; // netCDF group path to model data
-
- final NshmpNetcdfCoordinates nshmpCoords;
-
- private Logger logger;
-
- /**
- * @param nshmpNetcdfFilename
- * @param nshmpBaseGroup
- * @throws IOException
- */
- public NshmpNetcdfReader(String nshmpNetcdfFilename) throws IOException {
-
- this.nshmpNetcdfFilename = nshmpNetcdfFilename;
-
- // NSHMP 2018a data netCDF group path to data
- nshmpBaseGroup = "/CONUS/2018/v1.1";
-
- // Set netcdf-java logging level
- // WARNING is default
- // SEVERE disables warnings
- logger = Logger.getLogger("ucar");
- logger.setLevel(Level.SEVERE);
-
- // Open netCDF file
- // TODO: switch to try-with-resource block (solve errors when initializing
- // final fields)
- NetcdfDataset ncd = NetcdfDataset.openDataset(this.nshmpNetcdfFilename);
-
- // TODO: handle metadata from netCDF file
- // get root group and attributes, and attributes of subgroups
-
- // TODO: get netCDF dimensions
-
- // get target data group
- Group targetGroup = ncd.findGroup(nshmpBaseGroup);
-
- // get coordinate variables
- Variable vVs30 = targetGroup.findVariable("vs30");
- Variable vImts = targetGroup.findVariable("Imt");
- Variable vlats = targetGroup.findVariable("lat");
- Variable vlons = targetGroup.findVariable("lon");
- Variable vImls = targetGroup.findVariable("Imls");
- // TODO: hazard map variable will be renamed in netCDF file
- Variable vhazardLevels = targetGroup.findVariable("Poe");
- // or, for shared variables or dimensions:
- // Variable vLat = targetGroup.findVariableOrInParent("lat");
-
- // store coordinate variables in nshmpCoords
- nshmpCoords = new NshmpNetcdfCoordinates(vVs30, vImts, vlats, vlons, vImls, vhazardLevels);
-
- // Close NetcdfDataset
- // TODO: avoid this with try-with-resource block
- if (null != ncd) try {
- ncd.close();
- } catch (IOException ioe) {
- logger.log(Level.SEVERE, "Failed trying to close " + this.nshmpNetcdfFilename, ioe);
- }
-
- }
-
- /*
- * Return the full set of hazard curves at the bounding grid points and the
- * target Location. Intended for validation uses.
- */
- public static List<Map<SiteClass, Map<Imt, XySequence>>> getBoundingHazards(
- NshmpNetcdfReader nshmpNCRD, Location site) {
- BoundingHazards boundingHazards = new BoundingHazards(nshmpNCRD, site);
-
- return boundingHazards.boundingHazardsList;
- }
-
- /*
- * Return a Map<SiteClass, Map<Imt, XySequence>> with hazard curves at the
- * specified Location.
- */
- public static Map<SiteClass, Map<Imt, XySequence>> getHazardAtSite(NshmpNetcdfReader nshmpNCRD,
- Location site) {
- List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazardsList =
- getBoundingHazards(nshmpNCRD, site);
- return boundingHazardsList.get(4);
- }
-
- /*
- * Return a single hazard curve for the specified Location, SiteClass, and
- * Imt.
- */
- public static XySequence getHazardCurve(NshmpNetcdfReader nshmpNCRD, Location site, SiteClass sc,
- Imt imt) {
- Map<SiteClass, Map<Imt, XySequence>> hazardCurves = getHazardAtSite(nshmpNCRD, site);
- return hazardCurves.get(sc).get(imt);
- }
-
- /*
- * Container for dimensions and coordinate variables of nshmp netCDF file to
- * facilitate indexing operations prior to data retrieval.
- */
- private class NshmpNetcdfCoordinates {
-
- final int[] vs30s;
- final List<SiteClass> siteClassEnums;
-
- final double[] imtVals;
- final List<Imt> imts;
-
- final double[] lats;
- final double[] lons;
-
- final Map<SiteClass, Map<Imt, double[]>> imlMap;
- final double[] hazardLevels;
-
- final double maxLon;
- final double minLon;
- final double maxLat;
- final double minLat;
- final double deltaLon;
- final double deltaLat;
-
- final int nVs;
- final int nImt;
- final int nLat;
- final int nLon;
- final int nIml;
- final int nHazardLevel;
-
- /**
- * @param vVs30s
- * @param vImts
- * @param vlats
- * @param vlons
- * @param vImls
- * @param vHazardLevels
- * @throws IOException
- */
- public NshmpNetcdfCoordinates(
- Variable vVs30s,
- Variable vImts,
- Variable vlats,
- Variable vlons,
- Variable vImls,
- Variable vHazardLevels) throws IOException {
- // This bypasses the netCDF dimensions, but since we know what the
- // variables and their dimensions should be, this is OK(???)
-
- // read coordinate variables into java arrays
- vs30s = (int[]) vVs30s.read().get1DJavaArray(DataType.INT);
- imtVals = (double[]) vImts.read().get1DJavaArray(DataType.DOUBLE);;
- lats = (double[]) vlats.read().get1DJavaArray(DataType.DOUBLE);;
- lons = (double[]) vlons.read().get1DJavaArray(DataType.DOUBLE);;
- hazardLevels = (double[]) vHazardLevels.read().get1DJavaArray(DataType.DOUBLE);
-
- // assumes lons and lats are sorted ascending
- maxLon = lons[lons.length - 1];
- minLon = lons[0];
- maxLat = lats[lats.length - 1];
- minLat = lats[0];
-
- // get lat/lon spacing --> this is currently a required attribute but is
- // not used
- deltaLon = (double) vlons.findAttribute("spacing").getNumericValue();
- deltaLat = (double) vlats.findAttribute("spacing").getNumericValue();
-
- nVs = vs30s.length;
- nImt = imtVals.length;
- nLat = lats.length;
- nLon = lons.length;
- nHazardLevel = hazardLevels.length;
-
- // vImls has dimensions (vs30, Imt, Iml)
- // alternatively get nIml from Dimension Iml
- nIml = vImls.getDimension(2).getLength(); // Imls.length;
-
- // convert Vs30<int> to SiteClass Enum
- // TODO: site class should be stored as enum in netCDF file rather than
- // vs30 value, so this will break:
- siteClassEnums = Arrays.stream(vs30s)
- .mapToObj(x -> mapVsToEnum2018(x))
- .collect(Collectors.toUnmodifiableList());
-
- // convert Imt<double> to Imt Enum
- // TODO: Imt should be stored in netCDF as enum, so this will break:
- imts = Arrays.stream(imtVals)
- .mapToObj(x -> (x < 0.009) ? Imt.PGA : Imt.fromPeriod(x))
- .collect(Collectors.toUnmodifiableList());
-
- // Extract Imls to Map
- // vImls is a 3-dimensional array with dims (vs30, Imt, Iml)
- // Array ImlsArray = vImls.read();
- // int[] imlDims = ImlsArray.getShape();
-
- imlMap = mapImls(vImls);
- }
-
- /*
- * validate target coordinates
- */
- void checkCoords(double lon, double lat) {
- checkArgument(
- (lon >= minLon && lon <= maxLon && lat >= minLat && lat <= maxLat),
- String.format("Target coordinates out of range, %.3f <= lon < %.3f, %.3f <= lat <= %.3f",
- minLon, maxLat, minLat, maxLat));
- }
-
- void checkCoords(Location targetSite) {
- checkCoords(targetSite.longitude, targetSite.latitude);
- }
-
- /*
- * convert 3D Iml variable (dimensions vs30, Imt, Iml) to map of Imls by
- * SiteClass and Imt
- *
- * TODO: use MultiMap or SetMultiMap (etc.) to store unique IML sets? Could
- * then also initialize the underlying XySequence objects for reading in the
- * hazard curves...
- */
- private Map<SiteClass, Map<Imt, double[]>> mapImls(Variable vImls) {
-
- // initialize
- EnumMap<SiteClass, Map<Imt, double[]>> vsImtImlMap = Maps.newEnumMap(SiteClass.class);
-
- for (int i = 0; i < siteClassEnums.size(); i++) {
- SiteClass sc = siteClassEnums.get(i);
-
- Map<Imt, double[]> imtImlMap = Maps.newEnumMap(Imt.class);
- for (int j = 0; j < imts.size(); j++) {
- Imt imt = imts.get(j);
-
- // set origin and shape of double[] Imls to read
- int[] origin = new int[] { i, j, 0 };
- int[] shape = new int[] { 1, 1, nIml };
-
- try {
- imtImlMap.put(imt,
- (double[]) vImls.read(origin, shape).reduce().get1DJavaArray(DataType.DOUBLE));
- } catch (IOException | InvalidRangeException e) {
- System.err.println("Failed read attempt for vImls with origin: " +
- Arrays.toString(origin) + ", shape: " + Arrays.toString(shape));
- e.printStackTrace();
- }
- }
-
- vsImtImlMap.put(sc, imtImlMap);
- }
-
- return Maps.immutableEnumMap(vsImtImlMap);
- }
-
- private Location getLoc(int idxLon, int idxLat) {
- return Location.create(lats[idxLat], lons[idxLon]);
- }
- }
-
- /*
- * Container for gridded hazard curves at four closest grid points to target
- */
- static class BoundingHazards {
-
- final NshmpNetcdfReader ncrd;
- final Location site;
-
- private int idxLonLL;
- private int idxLatLL;
- private List<Location> boundingSites;
-
- final double fracLon;
- final double fracLat;
-
- final List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazardsList;
- // final Map<SiteClass, Map<Imt, double[]>> targetHazards;
-
- BoundingHazards(NshmpNetcdfReader ncrd, Location site) {
- this.ncrd = ncrd;
- this.site = site;
-
- // validate target coordinates
- ncrd.nshmpCoords.checkCoords(this.site);
-
- // Compute indices for lower left bounding grid point
- idxLonLL = getIdxLTEQ(ncrd.nshmpCoords.lons, this.site.longitude);
- idxLatLL = getIdxLTEQ(ncrd.nshmpCoords.lats, this.site.latitude);
-
- // Get bounding grid points: LL, UL, UR, LR
- boundingSites = new ArrayList<Location>(4);
- boundingSites.add(ncrd.nshmpCoords.getLoc(idxLonLL(), idxLatLL()));
- boundingSites.add(ncrd.nshmpCoords.getLoc(idxLonUL(), idxLatUL()));
- boundingSites.add(ncrd.nshmpCoords.getLoc(idxLonUR(), idxLatUR()));
- boundingSites.add(ncrd.nshmpCoords.getLoc(idxLonLR(), idxLatLR()));
-
- // Compute grid fractions to target point
- fracLon = calcGridFrac(ncrd.nshmpCoords.lons, idxLonLL(), site.longitude);
- fracLat = calcGridFrac(ncrd.nshmpCoords.lats, idxLatLL(), site.latitude);
-
- boundingHazardsList = this.ncrd.extractHazardsAt(idxLonLL, idxLatLL);
-
- // targetHazards = calcTargetHazards();
- boundingHazardsList.add(calcTargetHazards());
-
- }
-
- private Map<SiteClass, Map<Imt, XySequence>> calcTargetHazards() {
-
- // interpolate from LL to UL point
- Map<SiteClass, Map<Imt, XySequence>> westTarget =
- getTargetData(boundingHazardsList.get(0), boundingHazardsList.get(1), fracLat);
- // interpolate from LR to UR point
- Map<SiteClass, Map<Imt, XySequence>> eastTarget =
- getTargetData(boundingHazardsList.get(3), boundingHazardsList.get(2), fracLat);
-
- // do E-W interpolation
- return getTargetData(westTarget, eastTarget, fracLon);
-
- }
-
- void dumpCoords() {
- System.err.printf("%14s %14s %14s %14s %14s %14s %s\n", "Point", "LL", "UL", "UR", "LR",
- "Target", "Frac");
- System.err.printf("%14s %14.5f %14.5f %14.5f %14.5f %14.5f %.8f\n", "Lat:",
- boundingSites.get(0).latitude,
- boundingSites.get(1).latitude,
- boundingSites.get(2).latitude,
- boundingSites.get(3).latitude,
- site.latitude,
- fracLat);
- System.err.printf("%14s %14.5f %14.5f %14.5f %14.5f %14.5f %.8f\n", "Lon:",
- boundingSites.get(0).longitude,
- boundingSites.get(1).longitude,
- boundingSites.get(2).longitude,
- boundingSites.get(3).longitude,
- site.longitude,
- fracLon);
-
- }
-
- void dumpBoundingPlusTarget(SiteClass sc, Imt imt) {
- System.err.printf("Results for:\n");
- System.err.printf(" SiteClass: %2s, %d m/s\n", sc.toString(), mapSiteClassToVs30(sc));
- System.err.printf(" IMT: %s, %s\n", imt.name(), imt.period());
-
- dumpCoords();
- System.err.printf("%14s\n", "Iml");
- for (int i = 0; i < ncrd.nshmpCoords.nIml; i++) {
- // System.err.printf("%14d %14d %14d %14d %14d %14d\n",
- // (long) ncrd.nshmpCoords.imlMap.get(sc).get(imt)[i],
- // (long) boundingHazardsList.get(0).get(sc).get(imt).y(i),
- // (long) boundingHazardsList.get(1).get(sc).get(imt).y(i),
- // (long) boundingHazardsList.get(2).get(sc).get(imt).y(i),
- // (long) boundingHazardsList.get(3).get(sc).get(imt).y(i),
- // (long) boundingHazardsList.get(4).get(sc).get(imt).y(i));
- System.err.printf("%14.5f %14.6e %14.6e %14.6e %14.6e %14.6e\n",
- ncrd.nshmpCoords.imlMap.get(sc).get(imt)[i],
- boundingHazardsList.get(0).get(sc).get(imt).y(i),
- boundingHazardsList.get(1).get(sc).get(imt).y(i),
- boundingHazardsList.get(2).get(sc).get(imt).y(i),
- boundingHazardsList.get(3).get(sc).get(imt).y(i),
- boundingHazardsList.get(4).get(sc).get(imt).y(i));
- }
-
- }
-
- void dumpBoundingPlusTarget(int idxSc, int idxImt) {
- dumpBoundingPlusTarget(
- ncrd.nshmpCoords.siteClassEnums.get(idxSc),
- ncrd.nshmpCoords.imts.get(idxImt));
- }
-
- int idxLonLL() {
- return idxLonLL;
- }
-
- int idxLatLL() {
- return idxLatLL;
- }
-
- int idxLonUL() {
- return idxLonLL;
- }
-
- int idxLatUL() {
- return idxLatLL + 1;
- }
-
- int idxLonUR() {
- return idxLonLL + 1;
- }
-
- int idxLatUR() {
- return idxLatLL + 1;
- }
-
- int idxLonLR() {
- return idxLonLL + 1;
- }
-
- int idxLatLR() {
- return idxLatLL;
- }
-
- }
-
- /*
- * Return hazard curves for four closest grid points as a List of:
- * Map<SiteClass, Map<Imt, XySequence(iml,hazard)>>
- *
- * List order is clockwise from lower left corner: LL, UL, UR, LR, [T] with an
- * empty slot for the interpolated target hazards
- */
- private List<Map<SiteClass, Map<Imt, XySequence>>> extractHazardsAt(int idxLonLL, int idxLatLL) {
-
- // initialize list
- List<Map<SiteClass, Map<Imt, XySequence>>> boundingHazardMaps =
- new ArrayList<Map<SiteClass, Map<Imt, XySequence>>>(5);
-
- try (NetcdfDataset ncd = NetcdfDataset.openDataset(nshmpNetcdfFilename)) {
-
- Group targetGroup = ncd.findGroup(nshmpBaseGroup);
-
- // TODO: rename variable in netCDF
- Variable vHazards = targetGroup.findVariable("AEPs");
-
- // set up origin and shape arrays for reading hazard curves at four
- // bounding
- // (lon,lat) grid points
- int[] origin = new int[] { 0, 0, idxLatLL, idxLonLL, 0 };
- int[] shape = new int[] { nshmpCoords.nVs, nshmpCoords.nImt, 2, 2, nshmpCoords.nIml };
-
- // read data into Array for bounding grid points
- // System.err.println("reading hazards for four points");
- Array aHazards = vHazards.read(origin, shape);
-
- // Array aHazards now has shape [nVs,nImt,2,2,nIml]
- // ...so origin will now be [0,0,0,0,0] for LL grid point
- // ...and shape of requested array is [nVs,nImt,1,1,nIml]
-
- // System.err.println(aHazards.shapeToString());
- // Array aHazards is size 2 on the lon and lat dimensions
- shape[2] = 1;
- shape[3] = 1;
-
- // TODO: XySequence instead of double[]
-
- // Put data into hazard map for LL point
- origin[2] = 0; // lat index position
- origin[3] = 0; // lon index position
- // System.err.println("getting LL hazards");
- boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
- // Put data into hazard map for UL point
- origin[2] = 1; // lat index position
- origin[3] = 0; // lon index position
- // System.err.println("getting UL hazards");
- boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
- // Put data into hazard map for UR point
- origin[2] = 1; // lat index position
- origin[3] = 1; // lon index position
- // System.err.println("getting UR hazards");
- boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
- // Put data into hazard map for LR point
- origin[2] = 0; // lat index position
- origin[3] = 1; // lon index position
- // System.err.println("getting LR hazards");
- boundingHazardMaps.add(mapHazardsFromArray(aHazards.section(origin, shape)));
-
- // // read data into Array (using Variable.read() for each hazard curve)
- // THIS IS SLOW
- // // TODO: XySequence instead of double[]
- // // TODO: read data from variable ONCE!!! convert Array to maps
- // // initialize list
- // List<Map<SiteClass, Map<Imt, double[]>>> boundingHazardMaps =
- // new ArrayList<Map<SiteClass, Map<Imt, double[]>>>(4);
- // // Read data into hazard map for LL point
- // System.err.println("getting LL hazards");
- // boundingHazardMaps.add(mapHazards(vHazards, idxLonLL, idxLatLL));
- // // Read data into hazard map for UL point
- // System.err.println("getting UL hazards");
- // boundingHazardMaps.add(mapHazards(vHazards, idxLonLL, idxLatLL + 1));
- // // Read data into hazard map for UR point
- // System.err.println("getting UR hazards");
- // boundingHazardMaps.add(mapHazards(vHazards, idxLonLL + 1, idxLatLL +
- // 1));
- // // Read data into hazard map for LR point
- // System.err.println("getting LR hazards");
- // boundingHazardMaps.add(mapHazards(vHazards, idxLonLL + 1, idxLatLL));
-
- } catch (IOException e) {
- // TODO Auto-generated catch block
- e.printStackTrace();
- } catch (InvalidRangeException e) {
- // TODO Auto-generated catch block
- e.printStackTrace();
- }
-
- // interpolate hazards to target point
- return boundingHazardMaps;
-
- }
-
- /*
- * TODO: attempt to calculate index
- *
- * possibly as NshmpNetcdfCoordinates.calcIdxLTEQ()
- */
-
- /*
- * Read hazard curves from netCDF variable into map of hazards by SiteClass
- * and Imt
- *
- * TODO: if target is on a grid point (or on a grid lat or lon), no need to
- * read 4 bounding points ?
- */
- private Map<SiteClass, Map<Imt, XySequence>> mapHazardsFromArray(Array aHazards) {
- // aHazards has had its lat and lon dimensions reduced (so it now is rank 3)
- // aHazards[SiteClass, Imt, Iml]
-
- // initialize
- EnumMap<SiteClass, Map<Imt, XySequence>> vsImtHazardMap = Maps.newEnumMap(SiteClass.class);
-
- for (int iSc = 0; iSc < nshmpCoords.siteClassEnums.size(); iSc++) {
- SiteClass sc = nshmpCoords.siteClassEnums.get(iSc);
-
- Map<Imt, XySequence> imtHazardMap = Maps.newEnumMap(Imt.class);
- for (int iImt = 0; iImt < nshmpCoords.imts.size(); iImt++) {
- Imt imt = nshmpCoords.imts.get(iImt);
-
- // set origin and shape of double[] hazards to read
- int[] origin = new int[] { iSc, iImt, 0 };
- int[] shape = new int[] { 1, 1, nshmpCoords.nIml };
-
- // try block needed for sectioning Array
- try {
- imtHazardMap.put(imt,
- XySequence.create(
- nshmpCoords.imlMap.get(sc).get(imt),
- (double[]) aHazards.section(origin, shape).reduce()
- .get1DJavaArray(DataType.DOUBLE)));
- } catch (InvalidRangeException e) {
- // TODO Auto-generated catch block
- e.printStackTrace();
- }
-
- }
-
- vsImtHazardMap.put(sc, imtHazardMap);
- }
-
- return Maps.immutableEnumMap(vsImtHazardMap);
- }
-
- /*
- * Read hazard curves from netCDF variable into map of hazards by SiteClass
- * and Imt
- *
- * TODO: if target is on a grid point (or on a grid lat or lon), no need to
- * read 4 bounding points
- */
- private Map<SiteClass, Map<Imt, double[]>> mapHazards(Variable vHazards, int idxLon, int idxLat) {
-
- // initialize
- EnumMap<SiteClass, Map<Imt, double[]>> vsImtHazardMap = Maps.newEnumMap(SiteClass.class);
-
- for (int iSc = 0; iSc < nshmpCoords.siteClassEnums.size(); iSc++) {
- SiteClass sc = nshmpCoords.siteClassEnums.get(iSc);
-
- Map<Imt, double[]> imtHazardMap = Maps.newEnumMap(Imt.class);
- for (int iImt = 0; iImt < nshmpCoords.imts.size(); iImt++) {
- Imt imt = nshmpCoords.imts.get(iImt);
-
- // set origin and shape of double[] hazards to read
- int[] origin = new int[] { iSc, iImt, idxLat, idxLon, 0 };
- int[] shape = new int[] { 1, 1, 1, 1, nshmpCoords.nIml };
-
- try {
- imtHazardMap.put(imt,
- (double[]) vHazards.read(origin, shape).reduce().get1DJavaArray(DataType.DOUBLE));
- } catch (IOException | InvalidRangeException e) {
- System.err.println("Failed read attempt for vHazards with origin: " +
- Arrays.toString(origin) + ", shape: " + Arrays.toString(shape));
- e.printStackTrace();
- }
-
- }
-
- vsImtHazardMap.put(sc, imtHazardMap);
- }
-
- return Maps.immutableEnumMap(vsImtHazardMap);
- }
-
- /*
- * find index of first element in a (sorted ascending) that is less than or
- * equal to target value t. If target value is equal to the maximum value in a
- * (the last element), an index of a.length - 1 is returned.
- */
- private static int getIdxLTEQ(double[] a, double t) {
- // assumes array is in sorted order (a[i] < a[i+1])
- // make sure target is within the range of the array
- double tol = 0.000001;
- int n = a.length;
- if (t < a[0] || t > a[n - 1]) {
- // should never get here thanks to NshmpNetcdfCoordinates.checkCoords()
- throw new IllegalArgumentException(
- String.format("Target outside of valid range: %.4f <= t <= %4f", a[0], a[n - 1]));
- }
- // if (t == a[0]) {
- if (DoubleMath.fuzzyEquals(a[0], t, tol)) {
- return 0;
- }
- // if (t == a[n - 1]) {
- if (DoubleMath.fuzzyEquals(a[n - 1], t, tol)) {
- return n - 2; // return second to last index number
- }
-
- int idx = Arrays.binarySearch(a, t);
- if (idx < 0) {
- // "exact" match not found, so use index of first value less than target
- // this is insertion_point - 1
- // returned idx = -insertion_point - 1
- idx = -1 * (idx + 1) - 1;
- // array[idx] <= target
- }
- return idx;
- }
-
- /*
- * Calculate fractional distance from a1 to t, between a1 and a2
- */
- private static double calcFrac(double a1, double a2, double t) {
- double tol = 0.00001;
- double frac = 0.0;
- if (Math.abs(t - a1) < tol) {
- // target value == a1
- frac = 0.0;
- } else if (Math.abs(t - a2) < tol) {
- // target value == a2
- frac = 1.0;
- } else {
- // calculate fractional distance to t between a[i] and a[i+1]
- frac = (t - a1) / (a2 - a1);
- }
- return frac;
- }
-
- /*
- * Calculate fractional distance from a[i] to t, between a[i] and a[i+1]
- */
- private static double calcGridFrac(double[] a, int i, double t) {
- return calcFrac(a[i], a[i + 1], t);
- }
-
- /*
- * Get data for target point
- *
- * @param d1 data at first point (p1)
- *
- * @param d2 data at second point (p2)
- *
- * @param frac fractional distance between p1 and p2 to target point
- */
- private static Map<SiteClass, Map<Imt, XySequence>> getTargetData(
- Map<SiteClass, Map<Imt, XySequence>> d1,
- Map<SiteClass, Map<Imt, XySequence>> d2,
- double frac) {
- // do we need better checking here? or is it safe to assume that every
- // Map<SiteClass, Map<Imt,double[]>> passed in is consistent?
- if (d1.size() != d2.size()) {
- throw new IllegalArgumentException("Array size disagreement, cannot interpolate");
- }
- Map<SiteClass, Map<Imt, XySequence>> tMap;
- if (frac == 0.0) {
- // target is the same as d1
- tMap = d1;
- } else if (frac == 1.0) {
- // target is the same as d2
- tMap = d2;
- } else {
- tMap = linearInterpolate(d1, d2, frac);
- }
- return tMap;
- }
-
- /*
- * Linear interpolation of data values to a target point
- *
- * TODO: use XySequence.copyOf(), streams
- */
- private static Map<SiteClass, Map<Imt, XySequence>> linearInterpolate(
- Map<SiteClass, Map<Imt, XySequence>> v1,
- Map<SiteClass, Map<Imt, XySequence>> v2,
- double frac) {
- // return v1 * (1 - frac) + v2 * frac
- // DoubleData routines add and multiply in place, we don't want to modify
- // the original data here
- if (v1.size() != v2.size()) {
- throw new IllegalArgumentException("Array size disagreement, cannot interpolate");
- }
-
- Map<SiteClass, Map<Imt, XySequence>> tMap = Maps.newEnumMap(SiteClass.class);
-
- for (SiteClass sc : v1.keySet()) {
- Map<Imt, XySequence> imtHazards = Maps.newEnumMap(Imt.class);
-
- for (Imt imt : v1.get(sc).keySet()) {
- double[] v1Haz = v1.get(sc).get(imt).yValues().toArray();
- double[] v2Haz = v2.get(sc).get(imt).yValues().toArray();
- double[] t = new double[v1Haz.length];
- for (int i = 0; i < v1Haz.length; i++) {
- t[i] = v1Haz[i] * (1 - frac) + v2Haz[i] * frac;
- }
- imtHazards.put(imt, XySequence.create(v1.get(sc).get(imt).xValues().toArray(), t));
- }
- tMap.put(sc, imtHazards);
- }
- return tMap;
- }
-
- /*
- * check that sets of Iml are consistent for a given Imt across all SiteClass
- */
- private static boolean checkConsistentImls(Map<SiteClass, Map<Imt, double[]>> imlMap,
- List<SiteClass> siteClasses, List<Imt> imts) {
-
- boolean allEqual = true;
-
- // System.err.println("\n\n");
- for (Imt imt : imts) {
- // System.err.println("Checking Imt " + imt.toString());
-
- double[] firstImls = new double[1];
- SiteClass firstSc = null;
- for (SiteClass sc : siteClasses) {
- if (firstImls.length == 1) {
- firstImls = imlMap.get(sc).get(imt);
- firstSc = sc;
- continue;
- }
- boolean imlsDiffer = false;
- double[] theseImls = imlMap.get(sc).get(imt);
- for (int i = 0; i < firstImls.length; i++) {
- if (DoubleMath.fuzzyEquals(firstImls[i], theseImls[i], 1e-9)) {
- continue;
- }
- imlsDiffer = true;
- allEqual = false;
- break;
- }
- // if (imlsDiffer) {
- // // Imls differ for imt between lastSc and sc
- // System.err.println("Iml set differs, " + firstSc.toString() + " vs. "
- // + sc.toString());
- // System.err.println(" " + Arrays.toString(firstImls));
- // System.err.println(" " + Arrays.toString(theseImls));
- //
- // }
- }
- }
- return allEqual;
- }
-
- /**
- * This should be imported. For now, it is copied from
- * nshmp-haz-v2/src/gov.usgs.earthquake.nshmp.site.NehrpSiteClass.java
- *
- * Placeholder enum for likely move to Nehrp site class identifier instead of
- * Vs30.
- *
- * <p>These site class identifiers map to NEHRP site clases, but the intent is
- * that they can be used more generally for models in other parts of the world
- * where the GMMs are not necessarily parameterized in terms of vs30 to define
- * site response. For instance, NZ/JP site classes might use A, B, C, and D as
- * a proxy for the local I, II, III, and IV identifiers. In the U.S., models
- * will need to specify the Vs30 value that each site class corresponds to.
- * Although the values were consistent over prior models, now that multiple
- * site classes are supported (in 2018) across the entire U.S., there have
- * been changes proposed for balloting by the BSSC to make the Vs30
- * definitions of site classes consistent in how they are calculated.
- *
- * @author Peter Powers
- */
- enum SiteClass {
-
- /*
- * Notes on calculation of Vs30 for site class:
- *
- * Question: Why is it that the soil shear wave velocity shown in the
- * Unified Hazard Tool is not equal to the average of the values shown in
- * ASCE 7-10 table 20.3-1?
- *
- * For instance: 259 m/s (Site Class D), from the Unified Hazard Tool, is
- * not equal to (600 ft/s + 1200 ft/s)/2 * .3048 = 274 m/s
- *
- * Answer (Sanaz): we take the geometric mean: sqrt(1200*600)*0.3048 =
- * 258.6314 , which rounds to 259m/s.
- *
- *
- */
-
- /* OLD Vs30, NEW Vs30 */
-
- /* 2000 2000 */
- A,
-
- /* 1500 1500 (new) */
- AB,
-
- /* 1150 1080 */
- B,
-
- /* 760 760 */
- BC,
-
- /* 537 530 */
- C,
-
- /* 360 365 */
- CD,
-
- /* 259 260 */
- D,
-
- /* 180 185 */
- DE,
-
- /* 150 150 (new) */
- E,
-
- ;
- }
-
- static SiteClass mapVsToEnum2018(int vs) {
- // vs30 = 150, 185, 260, 365, 530, 760, 1080, 1500 ;
- switch (vs) {
- case 150:
- return SiteClass.E;
- case 185:
- return SiteClass.DE;
- case 260:
- return SiteClass.D;
- case 365:
- return SiteClass.CD;
- case 530:
- return SiteClass.C;
- case 760:
- return SiteClass.BC;
- case 1080:
- return SiteClass.B;
- case 1500:
- return SiteClass.AB;
- default:
- throw new IllegalArgumentException(
- String.format("Invalid Vs30, no matching 2018 SiteClass: %d", vs));
- }
- }
-
- static int mapSiteClassToVs30(SiteClass sc) {
- switch (sc) {
- case A:
- return 2000;
- case AB:
- return 1500;
- case B:
- return 1080;
- case BC:
- return 760;
- case C:
- return 530;
- case CD:
- return 365;
- case D:
- return 260;
- case DE:
- return 185;
- case E:
- return 150;
- default:
- throw new IllegalArgumentException(
- String.format("Invalid SiteClass %s", sc.toString()));
- }
- }
-}
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/netcdf/SiteClass.java b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/SiteClass.java
new file mode 100644
index 0000000000000000000000000000000000000000..89d22fef5415fa43f4004cb259ef8dbb0d5faf6e
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/netcdf/SiteClass.java
@@ -0,0 +1,98 @@
+package gov.usgs.earthquake.nshmp.netcdf;
+
+import java.util.Arrays;
+
+/**
+ * This should be imported. For now, it is copied from
+ * nshmp-haz-v2/src/gov.usgs.earthquake.nshmp.site.NehrpSiteClass.java
+ *
+ * Placeholder enum for likely move to Nehrp site class identifier instead of
+ * Vs30.
+ *
+ * <p>These site class identifiers map to NEHRP site clases, but the intent is
+ * that they can be used more generally for models in other parts of the world
+ * where the GMMs are not necessarily parameterized in terms of vs30 to define
+ * site response. For instance, NZ/JP site classes might use A, B, C, and D as a
+ * proxy for the local I, II, III, and IV identifiers. In the U.S., models will
+ * need to specify the Vs30 value that each site class corresponds to. Although
+ * the values were consistent over prior models, now that multiple site classes
+ * are supported (in 2018) across the entire U.S., there have been changes
+ * proposed for balloting by the BSSC to make the Vs30 definitions of site
+ * classes consistent in how they are calculated.
+ *
+ * @author Peter Powers
+ */
+public enum SiteClass {
+
+ /*
+ * Notes on calculation of Vs30 for site class:
+ *
+ * Question: Why is it that the soil shear wave velocity shown in the Unified
+ * Hazard Tool is not equal to the average of the values shown in ASCE 7-10
+ * table 20.3-1?
+ *
+ * For instance: 259 m/s (Site Class D), from the Unified Hazard Tool, is not
+ * equal to (600 ft/s + 1200 ft/s)/2 * .3048 = 274 m/s
+ *
+ * Answer (Sanaz): we take the geometric mean: sqrt(1200*600)*0.3048 =
+ * 258.6314 , which rounds to 259m/s.
+ *
+ *
+ */
+
+ /* OLD Vs30, NEW Vs30 */
+
+ /* 2000 2000 */
+ A(2000),
+
+ /* 1500 1500 (new) */
+ AB(1500),
+
+ /* 1150 1080 */
+ B(1080),
+
+ /* 760 760 */
+ BC(760),
+
+ /* 537 530 */
+ C(530),
+
+ /* 360 365 */
+ CD(365),
+
+ /* 259 260 */
+ D(260),
+
+ /* 180 185 */
+ DE(185),
+
+ /* 150 150 (new) */
+ E(150);
+
+ private final int vs30;
+
+ private SiteClass(int vs30) {
+ this.vs30 = vs30;
+ }
+
+ /**
+ * Returns the Vs30 associated with the {@code SiteClass}.
+ */
+ public int vs30() {
+ return vs30;
+ }
+
+ /**
+ * Returns a {@code SiteClass} associated with the vs30.
+ *
+ * @param vs30 The vs30 of the site class
+ */
+ public static SiteClass ofValue(int vs30) {
+ return Arrays.stream(values())
+ .filter(siteClass -> siteClass.vs30 == vs30)
+ .findFirst()
+ .orElseThrow(() -> new IllegalArgumentException(
+ "No matching site class to Vs30 [" + vs30 + "]"));
+ }
+
+}