//
// UsgsDemAdapter.java
//
/*
VisAD system for interactive analysis and visualization of numerical
data. Copyright (C) 1996 - 2017 Bill Hibbard, Curtis Rueden, Tom
Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and
Tommy Jasmin.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA
*/
package visad.data.gis;
import java.io.*;
import java.awt.geom.Rectangle2D;
import java.util.zip.GZIPInputStream;
import visad.*;
import visad.georef.*;
import visad.data.units.Parser;
import visad.util.ReflectedUniverse;
import java.util.Arrays;
import java.rmi.RemoteException;
/**
* UsgsDemAdapter converts the data in a USGS DEM file to a VisAD object.
* The specification for * the format can be found at
*
* http://rockyweb.cr.usgs.gov/nmpstds/acrodocs/dem/2DEM0198.PDF. Refer
* to this document for further explanation of the data elements of a DEM.
*
* To simplify things, class UsgsDem changes the format of a USGS 7.5 minute
* DEM to a simpler "row by column" raster form rather than the
* "row by column by subcolumn" format in the specification.
*
* @author Don Murray
* @author Jeffrey Albert Bergamini
* (http://srproj.lib.calpoly.edu/projects/csc/jbergam/reports/final.html)
* @version $Revision: 1.10 $ $Date: 2009-03-02 23:35:48 $
*/
public class UsgsDemAdapter {
/**************************************************************
* Fields related to java implementation
**************************************************************/
/** Southwest corner index constant */
private static final int SW = 0;
/** Northwest corner index constant */
private static final int NW = 1;
/** Northeast corner index constant */
private static final int NE = 2;
/** Southeast corner index constant */
private static final int SE = 3;
/** Geographic type reference system */
private static final int GEOGRAPHIC = 0;
/** UTM type reference system */
private static final int UTM = 1;
/** State type reference system */
private static final int STATE = 2;
/** Temporary storage string */
private String tempString;
/** Total number of points (vertices) in the DEM */
private int numPoints = 0;
/** A BufferedReader used to read from DEM files */
private BufferedReader in;
/** Marker to keep track of position within a record */
private int pos;
/** Record block size */
private static final int BLOCK_SIZE = 1024;
/** Temporary storage array */
private char[] array = new char[BLOCK_SIZE];
/**************************************************************
* Variables related to USGS DEM format
**************************************************************/
/** Measurement definition constant */
private static final int
RADIANS = 0,
FEET = 1,
METERS = 2,
ARCSECONDS = 3;
private static RealType EASTING = UTMCoordinate.EASTING;
private static RealType NORTHING = UTMCoordinate.NORTHING;
private static Unit[] measureUnits = new Unit[4];
static {
measureUnits[RADIANS] = CommonUnit.radian;
measureUnits[METERS] = CommonUnit.meter;
try {
measureUnits[FEET] = Parser.parse("foot");
measureUnits[ARCSECONDS] = Parser.parse("arcseconds");
} catch (Exception e) {}
}
private static int MISSING_VALUE = -32767;
/**
* File Name -- USGS DEM Logical Record Type A, Data Element 1
*
* The authorized digital cell name followed by a
* comma, space, and the two-character State
* designator(s) separated by hyphens. Abbreviations
* for other countries, such as Canada and Mexico,
* shall not be represented in the DEM header.
*/
private String fileName;
/**
* Free Format Text -- USGS DEM Logical Record Type A, Data Element 1
*
* Free format descriptor field, contains useful information
* related to digital process such as digitizing instrument,
* photo codes, slot widths, etc.
*/
private String freeFormatText;
/**
* Origin Code -- USGS DEM Logical Record Type A, Data Element 2
*
* 1=Autocorrelation RESAMPLE Simple bilinear
* 2=Manual profile GRIDEM Simple bilinear
* 3=DLG/hypsography CTOG 8-direction linear
* 4=Interpolation from photogrammetric system contours DCASS 4-direction linear
* 5=DLG/hypsography LINETRACE, LT4X Complex linear
* 6=DLG/hypsography CPS-3, ANUDEM, GRASS Complex polynomial
* 7=Electronic imaging (non-photogrametric), active or passive, sensor systems.
*
*/
private String originCode;
/**
* DEM Level Code -- USGS DEM Logical Record Type A, Data Element 3
*
* Code 1=DEM-1
* 2=DEM-2
* 3=DEM-3
* 4=DEM-4
*
*/
private int level;
/**
* Code defining ground planimetric reference system -- USGS DEM Logical Record Type A, Data Element 5
*
* Code 0=Geographic
* 1=UTM
* 2=State plane
* For codes 3-20, see Appendix 2-G.
* Code 0 represents the geographic
* (latitude/longitude) system for 30-minute,
* 1-degree and Alaska DEM's. Code 1 represents
* the current use of the UTM coordinate system for
* 7.5-minute DEM's
*
*/
private int referenceSystem;
/**
* Code defining zone in ground planimetric reference system.
* USGS DEM Logical Record Type A, Data Element 6
*/
private int zone;
/**
* Code defining unit of measure for ground planimetric coordinates.
* USGS DEM Logical Record Type A, Data Element 8
*
* 0 = radians
* 1 = feet
* 2 = meters
* 3 = arc-seconds
*/
private int groundMeasure;
/**
* Code defining unit of measure for elevation coordinates.
* USGS DEM Logical Record Type A, Data Element 9
*/
private int elevationMeasure;
/**
* A 4-element array containing the ground coordinates of the quadrangle
* boundary for the DEM.
*
* USGS DEM Logical Record Type A, Data Element 11
*/
private RealTuple[] corner = new RealTuple[4];
/** Minimum elevation */
private float minElevation;
/** Maximum elevation */
private float maxElevation;
/** Rotation angle elevation */
private float rotationAngle = 0;
/** Minimum easting */
private float minX = Float.POSITIVE_INFINITY;
/** Maximum easting */
private float maxX = Float.NEGATIVE_INFINITY;
/** Minimum northing */
private float minY = Float.POSITIVE_INFINITY;
/** Maximum northing */
private float maxY = Float.NEGATIVE_INFINITY;
/** Range of elevation */
private float elevationRange;
/** Range of easting */
private float xRange;
/** Range of northing */
private float yRange;
/**
* Elevation accuracy code
* USGS DEM Logical Record Type A, Data Element 14
*/
private int elevationAccuracy;
/**
* X (easting) spatial resolution
* USGS DEM Logical Record Type A, Data Element 15
*/
private float xResolution;
/**
* Y (northing) spatial resolution
* USGS DEM Logical Record Type A, Data Element 15
*/
private float yResolution;
/**
* Z (elevation) spatial resolution
* USGS DEM Logical Record Type A, Data Element 15
*/
private float zResolution;
/** Number of rows of profiles in the DEM */
// this is always 1 and superfluous
private int numRows;
/**
* Number of columns (in the northing dimension)
* of profiles in the DEM
*/
private int numColumns;
/**
* Maximum number of data points per column (easting dimension)
*/
private int maxRows;
/**
* float array of points in the form [xyz][numpoints]
*/
private float[][] rawCoords;
/**
* Vertical Datum code
*/
private int verticalDatum;
/**
* Horizontal Datum code
*/
private int horizontalDatum;
private Linear2DSet domainSet;
private FlatField data;
private boolean canDoGeotransform = true;
private ReflectedUniverse reflectUni;
/**
* Default constructor, everything remains null. Call load() to load a file.
*/
public UsgsDemAdapter () {
try {
reflectUni = new ReflectedUniverse();
reflectUni.exec("import ucar.visad.UTMCoordinateSystem");
reflectUni.exec("import geotransform.ellipsoids.Ellipsoid");
} catch (VisADException exc) { canDoGeotransform = false; }
}
/**
* Constructs a new UsgsDemAdapter object with data read from the given
* (native format, non-SDTS) USGS DEM file
*
* @param filename the name of the DEM file
*/
public UsgsDemAdapter (String filename) throws IOException, VisADException {
try {
reflectUni = new ReflectedUniverse();
reflectUni.exec("import ucar.visad.UTMCoordinateSystem");
reflectUni.exec("import geotransform.ellipsoids.Ellipsoid");
} catch (VisADException exc) { canDoGeotransform = false; }
if (filename != null) load(filename);
}
/**
* Reinitializes this UsgsDemAdapter object with data read from the given
* (non-SDTS) USGS DEM file.
*
* @param filename the name of the DEM file
* @throws IOException - If an I/O error occurs (invalid file)
*/
public void load(String filename) throws IOException, VisADException {
numPoints = 0;
data = null;
domainSet = null;
try {
//ungzip
if (filename.endsWith(".gz"))
in = new BufferedReader(
new InputStreamReader(
new GZIPInputStream(new FileInputStream(filename))));
else
in = new BufferedReader(new FileReader(filename));
} catch (Exception e) {
throw new IOException("Couldn't open file: " + filename);
}
processRecordTypeA();
for (int j=0; j= 0)
int m = parseInt(6); // 2 num elevations
int n = parseInt(6); // 2 n is always 1 -- hmm, useful
float Xgo = parseFloat(24); // x coord of starting point
float Ygo = parseFloat(24); // y coord of starting point
float localElevation = parseFloat(24); // elevation of local datum
float min = parseFloat(24); // min elevation for the profile
float max = parseFloat(24); // max elevation for the profile
if (m > maxRows) maxRows = m;
float costheta = (float) Math.cos(rotationAngle);
float sintheta = (float) Math.sin(rotationAngle);
/*
* Calculate profile coordinates
*/
int block = 0;
for (int j=0; j 0 && numToNextRecord < BLOCK_SIZE) {
skip(numToNextRecord);
}
// System.out.println("Size was " + pos);
// System.out.println("DONE WITH B");
}
/**
* Multi-line string representation
* @return An info string suitable for printing to the console.
*/
public String toString() {
StringBuffer s = new StringBuffer();
s.append(printElement("File name", fileName));
s.append(printElement("Free Format Text", freeFormatText));
s.append(printElement("Origin code", originCode));
s.append(printElement("Level code", level));
s.append(printElement("Reference system code", referenceSystem));
s.append(printElement("Reference system zone code", zone));
s.append(printElement("Ground measurement units", measureUnits[groundMeasure]));
s.append(printElement("Elev. measurement units", measureUnits[elevationMeasure]));
s.append(printElement("SW corner", corner[SW]));
s.append(printElement("NW corner", corner[NW]));
s.append(printElement("NE corner", corner[NE]));
s.append(printElement("SE corner", corner[SE]));
s.append(printElement("Min x", minX));
s.append(printElement("Max x", maxX));
s.append(printElement("Min y", minY));
s.append(printElement("Max y", maxY));
s.append(printElement("Min elevation", minElevation));
s.append(printElement("Max elevation", maxElevation));
s.append(printElement("Elevation accuracy code", elevationAccuracy));
s.append(printElement("X resolution", xResolution));
s.append(printElement("Y resolution", yResolution));
s.append(printElement("Z resolution", zResolution));
s.append(printElement("Columns of profiles", numColumns));
s.append(printElement("horizontal units", horizontalDatum));
s.append(printElement("vertical units", verticalDatum));
s.append(printElement("Max m in columns", maxRows));
s.append(printElement("Total data", numPoints));
return s.toString();
}
/**************************************************************
* Methods for gathering data elements from the file
**************************************************************/
/**
* Skips the specified number of bytes in the current file
*/
private void skip (int length) throws IOException {
try {
in.skip(length);
} catch (IOException e) {
throw new IOException("Could not skip " + length + " characters");
}
pos += length;
}
/**
* Reads in and returns a string of the specified length from the current
* position
*/
private String getString (int length) throws IOException {
if (length > 1024)
throw new IOException("Attempt to read more than 1024 bytes from file");
try {
in.read (array, 0, length);
} catch (IOException e) {
throw new IOException("Couldn't read string from file");
}
pos += length;
return (new String (array, 0 , length)).trim();
}
/**
* Reads in and returns a DEM-style integer from the current file position
* (taking up the specified number of bytes) as an int
*/
private int parseInt (int digits) throws IOException {
if (digits > 1024)
throw new IOException("Attempt to read more than 1024 bytes from file");
try {
in.read(array, 0, digits);
} catch (IOException e) {
throw new IOException("Couldn't read integer from file");
}
pos += digits;
// System.out.println("\"" + new String(array, 0, digits) + "\"");
return Integer.parseInt((new String(array, 0 , digits)).trim());
}
/**
* Reads in and returns a DEM-style integer of unknown length from the
* current position in the file as an int (necessary for non-standard DEMs)
*/
private int parseInt () throws IOException {
tempString = "";
char c;
try {
do {
c = (char) in.read();
pos++;
} while(Character.isWhitespace(c));
tempString += c;
do {
c = (char) in.read();
pos++;
tempString += c;
} while(!(Character.isWhitespace(c)));
} catch (IOException e) {
throw new IOException("Couldn't read integer from file");
}
tempString = tempString.substring(0,tempString.length()-1);
// System.out.println("tempString = " + tempString);
return Integer.parseInt(tempString);
}
/**
* Reads in and returns a DEM-style float from the current position in the
* file.
* @throws IOException If unable to read from file
*/
private float parseFloat (int digits) throws IOException {
//System.out.println("POS = " + pos);
if (digits > 1024)
throw new IOException("Attempt to read more than 1024 bytes from file");
try {
in.read(array, 0, digits);
} catch (IOException e) {
throw new IOException("Couldn't read integer from file");
}
pos += digits;
// System.out.println("\"" + new String(array, 0, digits) + "\"");
return translateReal((new String(array, 0 , digits)).trim());
}
/**
* Translates the String representation of a DEM-style real number
* into a float
*/
private float translateReal (String real) {
// System.out.println("\"" + real + "\"");
int plus, exp;
plus = real.indexOf('+');
if (plus > 0)
exp = Integer.parseInt(real.substring(plus+1,plus+3));
else
return Float.parseFloat(real);
if (real.charAt(0) != '-') {
tempString =
real.charAt(0) + real.substring(2,exp+2) + "." +
real.substring(exp+2,plus-1);
} else {
tempString = real.substring(0,2) + real.substring(3,exp+3) +
"." + real.substring(exp+2,plus-1);
}
return Float.parseFloat(tempString);
}
/**
* Loads a UTMCoordinate from the current position in the file
*/
private RealTuple parseCoordinate() throws IOException, VisADException {
if (referenceSystem == GEOGRAPHIC) {
// LatLonTuple is lat,lon order
float lon = parseFloat(24);
return (RealTuple) new LatLonTuple (parseFloat(24), lon);
} else if (referenceSystem == UTM) {
return (RealTuple)
new UTMCoordinate (parseFloat(24), parseFloat(24), zone);
} else {
throw new VisADException(
"Can't handle referenceSystem " + referenceSystem);
}
}
/** formatting stuff */
private String printElement(String title, Object value) {
return padRight(title,30) + ": " + value + "\n";
}
/** formatting stuff */
private String printElement(String title, int value) {
return printElement(title, new Integer(value));
}
/** formatting stuff */
private String printElement(String title, float value) {
return printElement(title, new Float(value));
}
/**
* This method takes any Object and using its String representation
* provided by its toString() method, pads it with blank characters
* on the right, to a specified length.
* @param obj Object to be padded
* @param i padding length
*/
public static String padRight(Object obj, int i) {
char filler = ' ';
String s = new String(obj.toString());
int j = s.length();
StringBuffer buf = new StringBuffer(s);
for(int k = 0; k < i - j; k++) buf.append(filler);
return buf.toString();
}
private Linear2DSet makeDomainSet() throws VisADException {
RealTupleType rtt = null;
Unit[] units = null;
switch (groundMeasure) {
case RADIANS:
units = new Unit[] {measureUnits[RADIANS], measureUnits[RADIANS]};
break;
case FEET:
units = new Unit[] {measureUnits[FEET], measureUnits[FEET]};
break;
case METERS:
units = new Unit[] {measureUnits[METERS], measureUnits[METERS]};
break;
case ARCSECONDS:
units = new Unit[] {measureUnits[ARCSECONDS], measureUnits[ARCSECONDS]};
break;
}
if (referenceSystem == UTM) {
// check to see if we can do uncompression
CoordinateSystem cs = null;
if (canDoGeotransform) {
try {
// the true says NH. Can we figure this out somehow?
switch (horizontalDatum) {
case 1:
reflectUni.setVar(
"ellipsoid", reflectUni.getVar("UTMCoordinateSystem.CC"));
break;
case 2:
reflectUni.setVar(
"ellipsoid", reflectUni.getVar("UTMCoordinateSystem.WD"));
break;
case 3:
reflectUni.setVar(
"ellipsoid", reflectUni.getVar("UTMCoordinateSystem.WE"));
break;
case 4:
reflectUni.setVar(
"ellipsoid", reflectUni.getVar("UTMCoordinateSystem.RF"));
break;
default:
reflectUni.setVar(
"ellipsoid", reflectUni.getVar("UTMCoordinateSystem.RF"));
}
reflectUni.setVar(
"bounds", new Rectangle2D.Double(minX, minY, maxX-minX, maxY-minY));
reflectUni.setVar("zone", zone);
reflectUni.setVar("north", true);
cs = (CoordinateSystem)
reflectUni.exec(
"cs = new UTMCoordinateSystem(ellipsoid, zone, north, bounds)");
cs = new CachingCoordinateSystem(cs);
} catch (VisADException ve) {
System.err.println(
"ucar.visad.UTMCoordinateSystem not found for UTM->lat/lon conversion" + ve);
cs = null;
}
rtt = new RealTupleType(EASTING, NORTHING, cs, null);
}
} else {
rtt = RealTupleType.SpatialEarth2DTuple;
}
domainSet =
(referenceSystem == UTM)
? new Linear2DSet(rtt,
//minX, maxX, ((int)(xRange/xResolution))+1,
minX, minX+(xResolution*(numColumns-1)), numColumns,
minY, maxY, ((int)(yRange/yResolution))+1,
(CoordinateSystem) null, /* in rtt if used */
units, (ErrorEstimate[]) null,
true /*cache*/)
/* GEOGRAPHIC */
: new LinearLatLonSet(rtt,
minX, minX+(xResolution*(numColumns-1)), numColumns,
minY, minY+(yResolution*(maxRows-1)), maxRows,
(CoordinateSystem) null, /* in rtt if used */
units, (ErrorEstimate[]) null,
true /*cache*/);
// System.out.println("domainSet = " + domainSet);
return domainSet;
}
private FlatField makeFlatField() throws VisADException {
makeDomainSet();
float[][] altitudes = new float[1][domainSet.getLength()];
Arrays.fill(altitudes[0], Float.NaN);
FunctionType ft =
new FunctionType(
((SetType)domainSet.getType()).getDomain(), RealType.Altitude);
FlatField ff = new FlatField(ft, domainSet,
(CoordinateSystem) null,
(Set[]) null,
new Unit[] {measureUnits[elevationMeasure]});
RealTuple coord = null;
int index = 0;
int[] indices =
domainSet.valueToIndex(new float[][] {rawCoords[0], rawCoords[1]});
float alt;
int numMissing = 0;
for (int i = 0; i < numPoints; i++) {
alt = rawCoords[2][i];
if (indices[i] < 0) {
numMissing++;
altitudes[0][indices[i]] = Float.NaN;
} else {
altitudes[0][indices[i]] = rawCoords[2][i];
}
}
// System.out.println("Out of " +numPoints+ " points, " + numMissing + " were missing");
try {
ff.setSamples(altitudes, false);
} catch (RemoteException re) {} // can't happen
data = ff;
return data;
}
/**
* Get the DEM as a VisAD data object
* @return a FlatField of type ((x,y) -> altitude) where x,y is either
* (Longitude, Latitude) or (EASTING, NORTHING)
*/
public FieldImpl getData() throws VisADException {
return (data == null) ? makeFlatField(): data;
}
/**
* Get the domain set for this DEM
* @return Gridded2DSet for domain
*/
public Gridded2DSet getDomain() throws VisADException {
return (domainSet == null) ? makeDomainSet() : domainSet;
}
public static void main(String[] args) throws Exception {
if (args.length == 0) {
System.out.println("Need to supply a filename");
System.exit(1);
}
UsgsDemAdapter uda = new UsgsDemAdapter(args[0]);
System.out.println(uda);
System.out.println(uda.getData().getDomainSet());
}
}
Ads help maintain this website.