ucar.nc2.dataset.conv.WRFConvention Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of cdm Show documentation
Show all versions of cdm Show documentation
The NetCDF-Java Library is a Java interface to NetCDF files,
as well as to many other types of scientific data formats.
The newest version!
/*
* Copyright 1998-2014 University Corporation for Atmospheric Research/Unidata
*
* Portions of this software were developed by the Unidata Program at the
* University Corporation for Atmospheric Research.
*
* Access and use of this software shall impose the following obligations
* and understandings on the user. The user is granted the right, without
* any fee or cost, to use, copy, modify, alter, enhance and distribute
* this software, and any derivative works thereof, and its supporting
* documentation for any purpose whatsoever, provided that this entire
* notice appears in all copies of the software, derivative works and
* supporting documentation. Further, UCAR requests that the user credit
* UCAR/Unidata in any publications that result from the use of this
* software or in any product that includes this software. The names UCAR
* and/or Unidata, however, may not be used in any advertising or publicity
* to endorse or promote any products or commercial entity unless specific
* written permission is obtained from UCAR/Unidata. The user also
* understands that UCAR/Unidata is not obligated to provide the user with
* any support, consulting, training or assistance of any kind with regard
* to the use, operation and performance of this software nor to provide
* the user with any updates, revisions, new versions or "bug fixes."
*
* THIS SOFTWARE IS PROVIDED BY UCAR/UNIDATA "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 UCAR/UNIDATA BE LIABLE FOR ANY SPECIAL,
* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
* FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
* NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
* WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package ucar.nc2.dataset.conv;
import ucar.ma2.*;
import ucar.nc2.*;
import ucar.nc2.constants.CDM;
import ucar.nc2.constants.CF;
import ucar.nc2.constants._Coordinate;
import ucar.nc2.constants.AxisType;
import ucar.nc2.time.*;
import ucar.nc2.units.SimpleUnit;
import ucar.nc2.util.CancelTask;
import ucar.nc2.dataset.*;
import ucar.nc2.dataset.transform.WRFEtaTransformBuilder;
import ucar.unidata.geoloc.*;
import ucar.unidata.geoloc.projection.*;
import ucar.unidata.util.StringUtil2;
import java.io.IOException;
import java.io.OutputStreamWriter;
import java.io.PrintWriter;
import java.util.*;
/**
* WRF netcdf output files.
*
* Note: Apparently WRF netcdf files before version 2 didnt output the projection origin, so
* we cant properly georeference them.
*
* This Convention currently only supports ARW output, identified as DYN_OPT=2 or GRIDTYPE=C
*
* @author caron
*/
public class WRFConvention extends CoordSysBuilder {
// static private java.text.SimpleDateFormat dateFormat;
// static {
// dateFormat = new java.text.SimpleDateFormat("yyyy-MM-dd_HH:mm:ss");
// dateFormat.setTimeZone(java.util.TimeZone.getTimeZone("GMT"));
//}
public static boolean isMine(NetcdfFile ncfile) {
if (null == ncfile.findDimension("south_north")) return false;
// ARW only
Attribute att = ncfile.findGlobalAttribute("DYN_OPT");
if (att != null) {
if (att.getNumericValue().intValue() != 2) return false;
} else {
att = ncfile.findGlobalAttribute("GRIDTYPE");
if (att == null) return false;
if (!att.getStringValue().equalsIgnoreCase("C") &&
!att.getStringValue().equalsIgnoreCase("E"))
return false;
}
att = ncfile.findGlobalAttribute("MAP_PROJ");
return att != null;
}
private double centerX = 0.0, centerY = 0.0;
private ProjectionCT projCT = null;
private boolean gridE = false;
public WRFConvention() {
this.conventionName = "WRF";
}
/* Implementation notes
There are 2 WRFs : NMM ("Non-hydrostatic Mesoscale Model" developed at NOAA/NCEP) and
1) NMM ("Non-hydrostatic Mesoscale Model" developed at NOAA/NCEP)
GRIDTYPE="E"
DYN_OPT = 4
This is a staggered grid that requires special processing.
2) ARW ("Advanced Research WRF", developed at MMM)
GRIDTYPE="C"
DYN_OPT = 2
DX, DY grid spaceing in meters (must be equal)
the Arakawa C staggered grid (see ARW 2.2 p 3-17)
the + are the "non-staggered" grid:
+ V + V + V +
U T U T U T U
+ V + V + V +
U T U T U T U
+ V + V + V +
U T U T U T U
+ V + V + V +
*/
/* ARW Users Guide p 3-19
7. MAP_PROJ_NAME: Character string specifying type of map projection. Valid entries are:
"polar" -> Polar stereographic
"lambert" -> Lambert conformal (secant and tangent)
"mercator" -> Mercator
8. MOAD_KNOWN_LAT/MOAD_KNOWN_LON (= CEN_LAT, CEN_LON):
Real latitude and longitude of the center point in the grid.
9. MOAD_STAND_LATS (= TRUE_LAT1/2):
2 real values for the "true" latitudes (where grid spacing is exact).
Must be between -90 and 90, and the values selected depend on projection:
Polar-stereographic: First value must be the latitude at which the grid
spacing is true. Most users will set this equal to their center latitude.
Second value must be +/-90. for NH/SH grids.
Lambert Conformal: Both values should have the same sign as the center
latitude. For a tangential lambert conformal, set both to the same value
(often equal to the center latitude). For a secant Lambert Conformal, they
may be set to different values.
Mercator: The first value should be set to the latitude you wish your grid
spacing to be true (often your center latitude). Second value is not used.
10. MOAD_STAND_LONS (=STAND_LON): This is one entry specifying the longitude in degrees East (-
180->180) that is parallel to the y-axis of your grid, (sometimes referred to as the
orientation of the grid). This should be set equal to the center longitude in most cases.
---------------------------------
version 3.1
http://www.mmm.ucar.edu/wrf/users/docs/user_guide_V3.1/users_guide_chap5.htm
The definition for map projection options:
map_proj = 1: Lambert Conformal
2: Polar Stereographic
3: Mercator
6: latitude and longitude (including global)
*/
public void augmentDataset(NetcdfDataset ds, CancelTask cancelTask) {
if (null != ds.findVariable("x")) return; // check if its already been done - aggregating enhanced datasets.
Attribute att = ds.findGlobalAttribute("GRIDTYPE");
gridE = att != null && att.getStringValue().equalsIgnoreCase("E");
// kludge in fixing the units
List vlist = ds.getVariables();
for (Variable v : vlist) {
att = v.findAttributeIgnoreCase(CDM.UNITS);
if (att != null) {
String units = att.getStringValue();
if (units != null)
v.addAttribute(new Attribute(CDM.UNITS, normalize(units))); // removes the old
}
}
// make projection transform
att = ds.findGlobalAttribute("MAP_PROJ");
int projType = att.getNumericValue().intValue();
boolean isLatLon = false;
if (projType == 203) {
/* centerX = centralLon;
centerY = centralLat;
ds.addCoordinateAxis( makeLonCoordAxis( ds, "longitude", ds.findDimension("west_east")));
ds.addCoordinateAxis( makeLatCoordAxis( ds, "latitude", ds.findDimension("south_north"))); */
Variable glat = ds.findVariable("GLAT");
if (glat == null) {
parseInfo.format("Projection type 203 - expected GLAT variable not found%n");
} else {
glat.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.Lat.toString()));
if (gridE) glat.addAttribute(new Attribute(_Coordinate.Stagger, CDM.ARAKAWA_E));
glat.setDimensions("south_north west_east");
glat.setCachedData(convertToDegrees(glat), false);
glat.addAttribute(new Attribute(CDM.UNITS, CDM.LAT_UNITS));
}
Variable glon = ds.findVariable("GLON");
if (glon == null) {
parseInfo.format("Projection type 203 - expected GLON variable not found%n");
} else {
glon.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.Lon.toString()));
if (gridE) glon.addAttribute(new Attribute(_Coordinate.Stagger, CDM.ARAKAWA_E));
glon.setDimensions("south_north west_east");
glon.setCachedData(convertToDegrees(glon), false);
glon.addAttribute(new Attribute(CDM.UNITS, CDM.LON_UNITS));
}
VariableDS v = new VariableDS(ds, null, null, "LatLonCoordSys", DataType.CHAR, "", null, null);
v.addAttribute(new Attribute(_Coordinate.Axes, "GLAT GLON Time"));
Array data = Array.factory(DataType.CHAR.getPrimitiveClassType(), new int[]{}, new char[]{' '});
v.setCachedData(data, true);
ds.addVariable(null, v);
Variable dataVar = ds.findVariable("LANDMASK");
dataVar.addAttribute(new Attribute(_Coordinate.Systems, "LatLonCoordSys"));
} else {
double lat1 = findAttributeDouble(ds, "TRUELAT1");
double lat2 = findAttributeDouble(ds, "TRUELAT2");
double centralLat = findAttributeDouble(ds, "CEN_LAT"); // center of grid
double centralLon = findAttributeDouble(ds, "CEN_LON"); // center of grid
double standardLon = findAttributeDouble(ds, "STAND_LON"); // true longitude
double standardLat = findAttributeDouble(ds, "MOAD_CEN_LAT");
ProjectionImpl proj = null;
switch (projType) {
case 0: // for diagnostic runs with no georeferencing
proj = new FlatEarth();
projCT = new ProjectionCT("flat_earth", "FGDC", proj);
// System.out.println(" using LC "+proj.paramsToString());
break;
case 1:
proj = new LambertConformal(standardLat, standardLon, lat1, lat2);
projCT = new ProjectionCT("Lambert", "FGDC", proj);
// System.out.println(" using LC "+proj.paramsToString());
break;
case 2:
// Thanks to Heiko Klein for figuring out WRF Stereographic
double lon0 = (Double.isNaN(standardLon)) ? centralLon : standardLon;
double lat0 = (Double.isNaN(centralLat)) ? lat2 : centralLat; // ?? 7/20/2010
double scaleFactor = (1 + Math.abs( Math.sin(Math.toRadians(lat1)))) / 2.; // R Schmunk 9/10/07
// proj = new Stereographic(lat2, lon0, scaleFactor);
proj = new Stereographic(lat0, lon0, scaleFactor);
projCT = new ProjectionCT("Stereographic", "FGDC", proj);
break;
case 3:
proj = new Mercator(standardLon, standardLat); // thanks to Robert Schmunk
projCT = new ProjectionCT("Mercator", "FGDC", proj);
// proj = new TransverseMercator(standardLat, standardLon, 1.0);
//projCT = new ProjectionCT("TransverseMercator", "FGDC", proj);
break;
case 6:
// version 3 "lat-lon", including global
// http://www.mmm.ucar.edu/wrf/users/workshops/WS2008/presentations/1-2.pdf
// use 2D XLAT, XLONG
isLatLon = true;
for (Variable v : vlist) {
if (v.getShortName().startsWith("XLAT")) {
v.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.Lat.toString()));
int[] shape = v.getShape();
if (v.getRank() == 3 && shape[0] == 1) { // remove time dependcies - MAJOR KLUDGE
List dims = v.getDimensions();
dims.remove(0);
v.setDimensions(dims);
}
} else if (v.getShortName().startsWith("XLONG")) {
v.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.Lon.toString()));
int[] shape = v.getShape();
if (v.getRank() == 3 && shape[0] == 1) { // remove time dependcies - MAJOR KLUDGE
List dims = v.getDimensions();
dims.remove(0);
v.setDimensions(dims);
}
} else if (v.getShortName().equals("T")) { // ANOTHER MAJOR KLUDGE to pick up 4D fields
v.addAttribute(new Attribute(_Coordinate.Axes, "Time XLAT XLONG z"));
} else if (v.getShortName().equals("U")) {
v.addAttribute(new Attribute(_Coordinate.Axes, "Time XLAT_U XLONG_U z"));
} else if (v.getShortName().equals("V")) {
v.addAttribute(new Attribute(_Coordinate.Axes, "Time XLAT_V XLONG_V z"));
} else if (v.getShortName().equals("W")) {
v.addAttribute(new Attribute(_Coordinate.Axes, "Time XLAT XLONG z_stag"));
}
}
break;
default:
parseInfo.format("ERROR: unknown projection type = %s%n", projType);
break;
}
if (proj != null) {
LatLonPointImpl lpt1 = new LatLonPointImpl(centralLat, centralLon); // center of the grid
ProjectionPoint ppt1 = proj.latLonToProj(lpt1, new ProjectionPointImpl());
centerX = ppt1.getX();
centerY = ppt1.getY();
if (debug) {
System.out.println("centerX=" + centerX);
System.out.println("centerY=" + centerY);
}
}
// make axes
if (!isLatLon) {
ds.addCoordinateAxis(makeXCoordAxis(ds, "x", ds.findDimension("west_east")));
ds.addCoordinateAxis(makeXCoordAxis(ds, "x_stag", ds.findDimension("west_east_stag")));
ds.addCoordinateAxis(makeYCoordAxis(ds, "y", ds.findDimension("south_north")));
ds.addCoordinateAxis(makeYCoordAxis(ds, "y_stag", ds.findDimension("south_north_stag")));
}
ds.addCoordinateAxis(makeZCoordAxis(ds, "z", ds.findDimension("bottom_top")));
ds.addCoordinateAxis(makeZCoordAxis(ds, "z_stag", ds.findDimension("bottom_top_stag")));
if (projCT != null) {
VariableDS v = makeCoordinateTransformVariable(ds, projCT);
v.addAttribute(new Attribute(_Coordinate.AxisTypes, "GeoX GeoY"));
if (gridE) v.addAttribute(new Attribute(_Coordinate.Stagger, CDM.ARAKAWA_E));
ds.addVariable(null, v);
}
}
// time coordinate variations
if (ds.findVariable("Time") == null) { // Can skip this if its already there, eg from NcML
CoordinateAxis taxis = makeTimeCoordAxis(ds, "Time", ds.findDimension("Time"));
if (taxis == null)
taxis = makeTimeCoordAxis(ds, "Time", ds.findDimension("Times"));
if (taxis != null)
ds.addCoordinateAxis(taxis);
}
ds.addCoordinateAxis(makeSoilDepthCoordAxis(ds, "ZS"));
ds.finish();
}
private Array convertToDegrees(Variable v) {
Array data;
try {
data = v.read();
data = data.reduce();
} catch (IOException ioe) {
throw new RuntimeException("data read failed on " + v.getFullName() + "=" + ioe.getMessage());
}
IndexIterator ii = data.getIndexIterator();
while (ii.hasNext()) {
ii.setDoubleCurrent(Math.toDegrees(ii.getDoubleNext()));
}
return data;
}
// pretty much WRF specific
private String normalize(String units) {
switch (units) {
case "fraction":
units = "";
break;
case "dimensionless":
units = "";
break;
case "NA":
units = "";
break;
case "-":
units = "";
break;
default:
units = StringUtil2.substitute(units, "**", "^");
units = StringUtil2.remove(units, '}');
units = StringUtil2.remove(units, '{');
break;
}
return units;
}
/////////////////////////////////////////////////////////////////////////
protected void makeCoordinateTransforms(NetcdfDataset ds) {
if (projCT != null) {
VarProcess vp = findVarProcess(projCT.getName(), null);
vp.isCoordinateTransform = true;
vp.ct = projCT;
}
super.makeCoordinateTransforms(ds);
}
protected AxisType getAxisType(NetcdfDataset ds, VariableEnhanced ve) {
Variable v = (Variable) ve;
String vname = v.getShortName();
if (vname.equalsIgnoreCase("x") || vname.equalsIgnoreCase("x_stag"))
return AxisType.GeoX;
if (vname.equalsIgnoreCase("lon"))
return AxisType.Lon;
if (vname.equalsIgnoreCase("y") || vname.equalsIgnoreCase("y_stag"))
return AxisType.GeoY;
if (vname.equalsIgnoreCase("lat"))
return AxisType.Lat;
if (vname.equalsIgnoreCase("z") || vname.equalsIgnoreCase("z_stag"))
return AxisType.GeoZ;
if (vname.equalsIgnoreCase("Z"))
return AxisType.Height;
if (vname.equalsIgnoreCase("time") || vname.equalsIgnoreCase("times"))
return AxisType.Time;
String unit = ve.getUnitsString();
if (unit != null) {
if (SimpleUnit.isCompatible("millibar", unit))
return AxisType.Pressure;
if (SimpleUnit.isCompatible("m", unit))
return AxisType.Height;
}
return null;
}
/**
* Does increasing values of Z go vertical up?
*
* @param v for thsi axis
* @return "up" if this is a Vertical (z) coordinate axis which goes up as coords get bigger,
* else return "down"
*/
public String getZisPositive(CoordinateAxis v) {
return "down"; //eta coords decrease upward
}
//////////////////////////////////////////////////////////////////////////////////////////////
private CoordinateAxis makeLonCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
double dx = 4 * findAttributeDouble(ds, "DX");
int nx = dim.getLength();
double startx = centerX - dx * (nx - 1) / 2;
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.DOUBLE, dim.getShortName(), "degrees_east", "synthesized longitude coordinate");
v.setValues(nx, startx, dx);
v.addAttribute(new Attribute(_Coordinate.AxisType, "Lon"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
return v;
}
private CoordinateAxis makeLatCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
double dy = findAttributeDouble(ds, "DY");
int ny = dim.getLength();
double starty = centerY - dy * (ny - 1) / 2;
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.DOUBLE, dim.getShortName(), "degrees_north", "synthesized latitude coordinate");
v.setValues( ny, starty, dy);
v.addAttribute(new Attribute(_Coordinate.AxisType, "Lat"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
return v;
}
private CoordinateAxis makeXCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
double dx = findAttributeDouble(ds, "DX") / 1000.0; // km ya just gotta know
int nx = dim.getLength();
double startx = centerX - dx * (nx - 1) / 2; // ya just gotta know
//System.out.println(" originX= "+originX+" startx= "+startx);
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.DOUBLE, dim.getShortName(), "km", "synthesized GeoX coordinate from DX attribute");
v.setValues(nx, startx, dx);
v.addAttribute(new Attribute(_Coordinate.AxisType, "GeoX"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
if (gridE) v.addAttribute(new Attribute(_Coordinate.Stagger, CDM.ARAKAWA_E));
return v;
}
private CoordinateAxis makeYCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
double dy = findAttributeDouble(ds, "DY") / 1000.0;
int ny = dim.getLength();
double starty = centerY - dy * (ny - 1) / 2; // - dy/2; // ya just gotta know
//System.out.println(" originY= "+originY+" starty= "+starty);
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.DOUBLE, dim.getShortName(), "km", "synthesized GeoY coordinate from DY attribute");
v.setValues(ny, starty, dy);
v.addAttribute(new Attribute(_Coordinate.AxisType, "GeoY"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
if (gridE) v.addAttribute(new Attribute(_Coordinate.Stagger, CDM.ARAKAWA_E));
return v;
}
private CoordinateAxis makeZCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
String fromWhere = axisName.endsWith("stag") ? "ZNW" : "ZNU";
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.DOUBLE, dim.getShortName(), "", "eta values from variable " + fromWhere);
v.addAttribute(new Attribute(CF.POSITIVE, CF.POSITIVE_DOWN)); // eta coordinate is 1.0 at bottom, 0 at top
v.addAttribute(new Attribute(_Coordinate.AxisType, "GeoZ"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
// create eta values from file variables: ZNU, ZNW
// But they are a function of time though the values are the same in the sample file
// NOTE: Use first time sample assuming all are the same!
Variable etaVar = ds.findVariable(fromWhere);
if (etaVar == null) return makeFakeCoordAxis(ds, axisName, dim);
int n = etaVar.getShape(1); //number of eta levels
int[] origin = new int[]{0, 0};
int[] shape = new int[]{1, n};
try {
Array array = etaVar.read(origin, shape);//read first time slice
ArrayDouble.D1 newArray = new ArrayDouble.D1(n);
IndexIterator it = array.getIndexIterator();
int count = 0;
while (it.hasNext()) {
double d = it.getDoubleNext();
newArray.set(count++, d);
}
v.setCachedData(newArray, true);
} catch (Exception e) {
e.printStackTrace();
}//ADD: error?
return v;
}
private CoordinateAxis makeFakeCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.SHORT, dim.getShortName(), "", "synthesized coordinate: only an index");
v.addAttribute(new Attribute(_Coordinate.AxisType, "GeoZ"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
v.setValues(dim.getLength(), 0, 1);
return v;
}
private CoordinateAxis makeTimeCoordAxis(NetcdfDataset ds, String axisName, Dimension dim) {
if (dim == null) return null;
int nt = dim.getLength();
Variable timeV = ds.findVariable("Times");
if (timeV == null) return null;
Array timeData;
try {
timeData = timeV.read();
} catch (IOException ioe) {
return null;
}
ArrayDouble.D1 values = new ArrayDouble.D1(nt);
int count = 0;
if (timeData instanceof ArrayChar) {
ArrayChar.StringIterator iter = ((ArrayChar) timeData).getStringIterator();
while (iter.hasNext()) {
String dateS = iter.next();
try {
CalendarDate cd = CalendarDateFormatter.isoStringToCalendarDate(null, dateS);
values.set(count++, (double) cd.getMillis() / 1000);
} catch (IllegalArgumentException e) {
parseInfo.format("ERROR: cant parse Time string = <%s> err= %s%n", dateS, e.getMessage());
// one more try
String startAtt = ds.findAttValueIgnoreCase(null, "START_DATE", null);
if ((nt == 1) && (null != startAtt)) {
try {
CalendarDate cd = CalendarDateFormatter.isoStringToCalendarDate(null, startAtt);
values.set(0, (double) cd.getMillis() / 1000);
} catch (IllegalArgumentException e2) {
parseInfo.format("ERROR: cant parse global attribute START_DATE = <%s> err=%s%n", startAtt, e2.getMessage());
}
}
}
}
} else {
IndexIterator iter = timeData.getIndexIterator();
while (iter.hasNext()) {
String dateS = (String) iter.next();
try {
CalendarDate cd = CalendarDateFormatter.isoStringToCalendarDate(null, dateS);
values.set(count++, (double) cd.getMillis() / 1000);
} catch (IllegalArgumentException e) {
parseInfo.format("ERROR: cant parse Time string = %s%n", dateS);
}
}
}
CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, DataType.DOUBLE, dim.getShortName(),
"secs since 1970-01-01 00:00:00", "synthesized time coordinate from Times(time)");
v.addAttribute(new Attribute(_Coordinate.AxisType, "Time"));
if (!axisName.equals(dim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, dim.getShortName()));
v.setCachedData(values, true);
return v;
}
private VariableDS makeSoilDepthCoordAxis(NetcdfDataset ds, String coordVarName) {
Variable coordVar = ds.findVariable(coordVarName);
if (null == coordVar)
return null;
Dimension soilDim = null;
List dims = coordVar.getDimensions();
for (Dimension d : dims) {
if (d.getShortName().startsWith("soil_layers"))
soilDim = d;
}
if (null == soilDim)
return null;
if (coordVar.getRank() == 1) {
coordVar.addAttribute(new Attribute(CF.POSITIVE, CF.POSITIVE_DOWN)); // soil depth gets larger as you go down
coordVar.addAttribute(new Attribute(_Coordinate.AxisType, "GeoZ"));
if (!coordVarName.equals(soilDim.getShortName()))
coordVar.addAttribute(new Attribute(_Coordinate.AliasForDimension, soilDim.getShortName()));
return (VariableDS) coordVar;
}
String units = ds.findAttValueIgnoreCase(coordVar, CDM.UNITS, "");
CoordinateAxis v = new CoordinateAxis1D(ds, null, "soilDepth", DataType.DOUBLE, soilDim.getShortName(), units, "soil depth");
v.addAttribute(new Attribute(CF.POSITIVE, CF.POSITIVE_DOWN)); // soil depth gets larger as you go down
v.addAttribute(new Attribute(_Coordinate.AxisType, "GeoZ"));
v.addAttribute(new Attribute(CDM.UNITS, CDM.UNITS));
if (!v.getShortName().equals(soilDim.getShortName()))
v.addAttribute(new Attribute(_Coordinate.AliasForDimension, soilDim.getShortName()));
//read first time slice
int n = coordVar.getShape(1);
int[] origin = new int[]{0, 0};
int[] shape = new int[]{1, n};
try {
Array array = coordVar.read(origin, shape);
ArrayDouble.D1 newArray = new ArrayDouble.D1(n);
IndexIterator it = array.getIndexIterator();
int count = 0;
while (it.hasNext()) {
double d = it.getDoubleNext();
newArray.set(count++, d);
}
v.setCachedData(newArray, true);
} catch (Exception e) {
e.printStackTrace();
}
return v;
}
private double findAttributeDouble(NetcdfDataset ds, String attname) {
Attribute att = ds.findGlobalAttributeIgnoreCase(attname);
if (att == null) return Double.NaN;
return att.getNumericValue().doubleValue();
}
/**
* Assign CoordinateTransform objects to Coordinate Systems.
*/
protected void assignCoordinateTransforms(NetcdfDataset ncDataset) {
super.assignCoordinateTransforms(ncDataset);
// any cs whose got a vertical coordinate with no units
List csys = ncDataset.getCoordinateSystems();
for (CoordinateSystem cs : csys) {
if (cs.getZaxis() != null) {
String units = cs.getZaxis().getUnitsString();
if ((units == null) || (units.trim().length() == 0)) {
VerticalCT vct = makeWRFEtaVerticalCoordinateTransform(ncDataset, cs);
if (vct != null)
cs.addCoordinateTransform(vct);
parseInfo.format("***Added WRFEta verticalCoordinateTransform to %s%n", cs.getName());
}
}
}
}
private VerticalCT makeWRFEtaVerticalCoordinateTransform(NetcdfDataset ds, CoordinateSystem cs) {
if ((null == ds.findVariable("PH")) || (null == ds.findVariable("PHB")) ||
(null == ds.findVariable("P")) || (null == ds.findVariable("PB")))
return null;
WRFEtaTransformBuilder builder = new WRFEtaTransformBuilder(cs);
return (VerticalCT) builder.makeCoordinateTransform(ds, null);
}
public static void main(String args[]) throws IOException, InvalidRangeException {
NetcdfFile ncd = NetcdfDataset.openFile("R:/testdata/wrf/WRFOU~C@", null);
Variable glat = ncd.findVariable("GLAT");
Array glatData = glat.read();
IndexIterator ii = glatData.getIndexIterator();
while (ii.hasNext()) {
ii.setDoubleCurrent(Math.toDegrees(ii.getDoubleNext()));
}
PrintWriter pw = new PrintWriter( new OutputStreamWriter(System.out, CDM.utf8Charset));
NCdumpW.printArray(glatData, "GLAT", pw, null);
Variable glon = ncd.findVariable("GLON");
Array glonData = glon.read();
ii = glonData.getIndexIterator();
while (ii.hasNext()) {
ii.setDoubleCurrent(Math.toDegrees(ii.getDoubleNext()));
}
NCdumpW.printArray(glonData, "GLON", pw, null);
Index index = glatData.getIndex();
Index index2 = glatData.getIndex();
int[] vshape = glatData.getShape();
int ny = vshape[1];
int nx = vshape[2];
ArrayDouble.D1 diff_y = (ArrayDouble.D1) Array.factory(DataType.DOUBLE, new int[]{ny});
ArrayDouble.D1 diff_x = (ArrayDouble.D1) Array.factory(DataType.DOUBLE, new int[]{nx});
for (int y = 0; y < ny - 1; y++) {
double val = glatData.getDouble(index.set(0, y, 0)) - glatData.getDouble(index2.set(0, y + 1, 0));
diff_y.set(y, val);
}
for (int x = 0; x < nx - 1; x++) {
double val = glatData.getDouble(index.set(0, 0, x)) - glatData.getDouble(index2.set(0, 0, x + 1));
diff_x.set(x, val);
}
NCdumpW.printArray(diff_y, "diff_y", pw, null);
NCdumpW.printArray(diff_x, "diff_x", pw, null);
ncd.close();
}
}
/*
NMM . output from gridgen. these are the input files to WRF i think
netcdf Q:/grid/netcdf/wrf/geo_nmm.d01.nc {
dimensions:
Time = UNLIMITED; // (1 currently)
DateStrLen = 19;
west_east = 136;
south_north = 309;
land_cat = 24;
soil_cat = 16;
month = 12;
variables:
char Times(Time=1, DateStrLen=19);
float XLAT_M(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "degrees latitude";
:description = "Latitude on mass grid";
:stagger = "M";
float XLONG_M(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "degrees longitude";
:description = "Longitude on mass grid";
:stagger = "M";
float XLAT_V(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "degrees latitude";
:description = "Latitude on velocity grid";
:stagger = "V";
float XLONG_V(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "degrees longitude";
:description = "Longitude on velocity grid";
:stagger = "V";
float E(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "-";
:description = "Coriolis E parameter";
:stagger = "M";
float F(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "-";
:description = "Coriolis F parameter";
:stagger = "M";
float LANDMASK(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "none";
:description = "Landmask : 1=land, 0=water";
:stagger = "M";
float LANDUSEF(Time=1, land_cat=24, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:units = "category";
:description = "24-category USGS landuse";
:stagger = "M";
float LU_INDEX(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "category";
:description = "Dominant category";
:stagger = "M";
float HGT_M(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "meters MSL";
:description = "Topography height";
:stagger = "M";
float HGT_V(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "meters MSL";
:description = "Topography height";
:stagger = "V";
float SOILTEMP(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "Kelvin";
:description = "Annual mean deep soil temperature";
:stagger = "M";
float SOILCTOP(Time=1, soil_cat=16, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:units = "category";
:description = "16-category top-layer soil type";
:stagger = "M";
float SOILCBOT(Time=1, soil_cat=16, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:units = "category";
:description = "16-category bottom-layer soil type";
:stagger = "M";
float ALBEDO12M(Time=1, month=12, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:units = "percent";
:description = "Monthly surface albedo";
:stagger = "M";
float GREENFRAC(Time=1, month=12, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XYZ";
:units = "fraction";
:description = "Monthly green fraction";
:stagger = "M";
float SNOALB(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "percent";
:description = "Maximum snow albedo";
:stagger = "M";
float SLOPECAT(Time=1, south_north=309, west_east=136);
:FieldType = 104; // int
:MemoryOrder = "XY ";
:units = "category";
:description = "Dominant category";
:stagger = "M";
:TITLE = "OUTPUT FROM GRIDGEN";
:SIMULATION_START_DATE = "0000-00-00_00:00:00";
:WEST-EAST_GRID_DIMENSION = 136; // int
:SOUTH-NORTH_GRID_DIMENSION = 309; // int
:BOTTOM-TOP_GRID_DIMENSION = 0; // int
:WEST-EAST_PATCH_START_UNSTAG = 1; // int
:WEST-EAST_PATCH_END_UNSTAG = 136; // int
:WEST-EAST_PATCH_START_STAG = 1; // int
:WEST-EAST_PATCH_END_STAG = 136; // int
:SOUTH-NORTH_PATCH_START_UNSTAG = 1; // int
:SOUTH-NORTH_PATCH_END_UNSTAG = 309; // int
:SOUTH-NORTH_PATCH_START_STAG = 1; // int
:SOUTH-NORTH_PATCH_END_STAG = 309; // int
:GRIDTYPE = "E";
:DX = 0.111143f; // float
:DY = 0.106776f; // float
:DYN_OPT = 4; // int
:CEN_LAT = 21.0f; // float
:CEN_LON = 80.0f; // float
:TRUELAT1 = 1.0E20f; // float
:TRUELAT2 = 1.0E20f; // float
:MOAD_CEN_LAT = 0.0f; // float
:STAND_LON = 1.0E20f; // float
:POLE_LAT = 90.0f; // float
:POLE_LON = 0.0f; // float
:corner_lats = 3.8832555f, 36.602543f, 36.602543f, 3.8832555f, 3.8931324f, 36.61482f, 36.59018f, 3.873307f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f; // float
:corner_lons = 65.589096f, 61.982986f, 98.01701f, 94.410904f, 65.69548f, 62.114895f, 98.14889f, 94.51727f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f; // float
:MAP_PROJ = 203; // int
:MMINLU = "USGS";
:ISWATER = 16; // int
:ISICE = 24; // int
:ISURBAN = 1; // int
:ISOILWATER = 14; // int
:grid_id = 1; // int
:parent_id = 1; // int
:i_parent_start = 0; // int
:j_parent_start = 0; // int
:i_parent_end = 136; // int
:j_parent_end = 309; // int
:parent_grid_ratio = 1; // int
}
*/ © 2015 - 2024 Weber Informatics LLC | Privacy Policy