ucar.nc2.dt.grid.CFGridWriter2 Maven / Gradle / Ivy
The newest version!
/*
* Copyright (c) 1998-2018 John Caron and University Corporation for Atmospheric Research/Unidata
* See LICENSE for license information.
*/
package ucar.nc2.dt.grid;
import ucar.ma2.Array;
import ucar.ma2.DataType;
import ucar.ma2.InvalidRangeException;
import ucar.ma2.Range;
import ucar.nc2.*;
import ucar.nc2.constants.*;
import ucar.nc2.dataset.*;
import ucar.nc2.dataset.transform.AbstractTransformBuilder;
import ucar.nc2.dt.GridCoordSystem;
import ucar.nc2.dt.GridDatatype;
import ucar.nc2.time.CalendarDate;
import ucar.nc2.time.CalendarDateRange;
import ucar.unidata.geoloc.*;
import ucar.unidata.geoloc.projection.LatLonProjection;
import ucar.unidata.util.Parameter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
/**
* Write CF compliant gridded data version 2.
*
* @deprecated Convert to CoverageCollection and use {@link ucar.nc2.ft2.coverage.writer.CFGridCoverageWriter}
*/
@Deprecated
public class CFGridWriter2 {
private static final org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(CFGridWriter2.class);
/**
* Compute the size of the file without writing
*
* @param gds the GridDataset
* @param gridList the list of variables to be written, or null for all
* @param llbb the lat/lon bounding box, or null for all
* @param projRect the projection bounding box, or null for all
* @param horizStride the x and y stride
* @param zRange the z stride
* @param dateRange date range, or null for all
* @param stride_time the time stride
* @param addLatLon add 2D lat/lon coordinates if needed
* @return total bytes written
*/
public static long makeSizeEstimate(ucar.nc2.dt.GridDataset gds, List gridList, LatLonRect llbb,
ProjectionRect projRect, int horizStride, Range zRange, CalendarDateRange dateRange, int stride_time,
boolean addLatLon) throws IOException, InvalidRangeException {
CFGridWriter2 writer2 = new CFGridWriter2();
return writer2.writeOrTestSize(gds, gridList, llbb, projRect, horizStride, zRange, dateRange, stride_time,
addLatLon, true, null);
}
/**
* Write a netcdf/CF file from a GridDataset
*
* @param gds the GridDataset
* @param gridList the list of variables to be written, or null for all
* @param llbb the lat/lon bounding box, or null for all
* @param projRect the projection bounding box, or null for all
* @param horizStride the x and y stride
* @param zRange the z stride
* @param dateRange date range, or null for all
* @param stride_time the time stride
* @param addLatLon add 2D lat/lon coordinates if needed
* @param writer this does the actual writing
* @return total bytes written
*/
public static long writeFile(ucar.nc2.dt.GridDataset gds, List gridList, LatLonRect llbb,
ProjectionRect projRect, int horizStride, Range zRange, CalendarDateRange dateRange, int stride_time,
boolean addLatLon, NetcdfFileWriter writer) throws IOException, InvalidRangeException {
CFGridWriter2 writer2 = new CFGridWriter2();
return writer2.writeOrTestSize(gds, gridList, llbb, projRect, horizStride, zRange, dateRange, stride_time,
addLatLon, false, writer);
}
/**
*
* @param gds the GridDataset
* @param gridList the list of variables to be written, or null for all
* @param llbb the lat/lon bounding box, or null for all
* @param projRect the projection bounding box, or null for all
* @param horizStride the x and y stride
* @param zRange the z stride
* @param dateRange date range, or null for all
* @param stride_time the time stride
* @param addLatLon add 2D lat/lon coordinates if needed
* @param testSizeOnly dont write, just return size
* @param writer this does the actual writing
* @return total bytes written
*/
private long writeOrTestSize(ucar.nc2.dt.GridDataset gds, List gridList, LatLonRect llbb,
ProjectionRect projRect, int horizStride, Range zRange, CalendarDateRange dateRange, int stride_time,
boolean addLatLon, boolean testSizeOnly, NetcdfFileWriter writer) throws IOException, InvalidRangeException {
NetcdfDataset ncd = (NetcdfDataset) gds.getNetcdfFile();
List varList = new ArrayList<>(); // could make these Sets
List varNameList = new ArrayList<>();
List axisList = new ArrayList<>();
if (gridList == null) { // want all of them
gridList = new ArrayList<>();
for (GridDatatype grid : gds.getGrids())
gridList.add(grid.getName());
}
LatLonRect resultBB = null;
long total_size = 0;
// add each desired Grid to the new file LOOK Might be better to group by coordSys
for (String gridName : gridList) {
if (varNameList.contains(gridName)) // already contains
continue;
GridDatatype gridOrg = gds.findGridDatatype(gridName);
if (gridOrg == null) {
log.debug("writeOrTestSize cant find grid {} - skipping%n", gridName);
continue;
}
GridCoordSystem gcsOrg = gridOrg.getCoordinateSystem();
// boolean isGlobal = gcsOrg.isGlobalLon();
varNameList.add(gridName);
// make subset as needed
Range timeRange = makeTimeRange(dateRange, gcsOrg.getTimeAxis1D(), stride_time);
Range zRangeUse = makeVerticalRange(zRange, gcsOrg.getVerticalAxis()); // LOOK no vert stride
GridDatatype gridWant;
List yxRanges = new ArrayList<>(2);
if (projRect != null) {
makeHorizRange(gcsOrg, projRect, horizStride, yxRanges);
gridWant = gridOrg.makeSubset(null, null, timeRange, zRangeUse, yxRanges.get(0), yxRanges.get(1));
} else {
if (llbb == null) {
yxRanges.add(null);
yxRanges.add(null);
} else {
yxRanges = gcsOrg.getRangesFromLatLonRect(llbb);
}
if ((null != timeRange) || (zRangeUse != null) || (llbb != null) || (horizStride > 1)) {
gridWant = gridOrg.makeSubset(timeRange, zRangeUse, llbb, 1, horizStride, horizStride);
} else
gridWant = gridOrg;
}
GridCoordSystem gcsWant = gridWant.getCoordinateSystem();
LatLonRect gridBB = gcsWant.getLatLonBoundingBox();
if (resultBB == null)
resultBB = gridBB;
else
resultBB.extend(gridBB);
Variable gridV = gridWant.getVariable(); // LOOK WTF ??
varList.add(gridV);
total_size += gridV.getSize() * gridV.getElementSize();
// add coordinate axes
addCoordinateAxis(gcsWant, varNameList, varList, axisList);
// add coordinate transform variables
addCoordinateTransform(gcsWant, ncd, varNameList, varList);
// Add Variables from the formula_terms
total_size += processTransformationVars(varList, varNameList, ncd, gds, gridWant, timeRange, zRangeUse,
yxRanges.get(0), yxRanges.get(1), horizStride, horizStride);
// optional lat/lon
if (addLatLon) {
Projection proj = gcsWant.getProjection();
if ((null != proj) && !(proj instanceof LatLonProjection)) {
total_size += addLatLon2D(ncd, varList, proj, gcsWant.getXHorizAxis(), gcsWant.getYHorizAxis());
addLatLon = false; // ??
}
}
}
if (testSizeOnly)
return total_size;
////////////////////////////////////////////////////////////////////
// start writing here
// check size is ok
boolean isLargeFile = isLargeFile(total_size);
writer.setLargeFile(isLargeFile);
if (resultBB != null)
addGlobalAttributes(writer, gds, resultBB);
// use fileWriter to copy the variables
FileWriter2 fileWriter = new FileWriter2(writer);
for (Variable v : varList)
fileWriter.addVariable(v);
addCFAnnotations(writer, gds, gridList, ncd, axisList, addLatLon);
writer.create();
// use fileWriter to copy the data
fileWriter.copyVarData(varList, null, null);
// updateGeospatialRanges(writer, llrect );
writer.close();
// this writes the data to the new file.
return total_size; // ok
}
private boolean isLargeFile(long total_size) {
boolean isLargeFile = false;
long maxSize = 2 * 1000 * 1000 * 1000; // LOOK why not use exact
if (total_size > maxSize) {
log.debug("Request size = {} Mbytes", total_size / 1000 / 1000);
isLargeFile = true;
}
return isLargeFile;
}
private void makeHorizRange(GridCoordSystem gcsOrg, ProjectionRect projRect, int horizStride, List yxRanges)
throws InvalidRangeException {
if (gcsOrg.getXHorizAxis().getRank() > 1 || gcsOrg.getYHorizAxis().getRank() > 1) {
throw new IllegalArgumentException("Coordinate systems with 2D horizontal axis are not supported");
}
CoordinateAxis1D xAxis = (CoordinateAxis1D) gcsOrg.getXHorizAxis();
double[] xCoords = xAxis.getCoordValues();
CoordinateAxis1D yAxis = (CoordinateAxis1D) gcsOrg.getYHorizAxis();
double[] yCoords = yAxis.getCoordValues();
ProjectionRect fullBB =
new ProjectionRect(xCoords[0], yCoords[0], xCoords[xCoords.length - 1], yCoords[yCoords.length - 1]);
if (!projRect.intersects(fullBB)) {
throw new InvalidRangeException("BBOX must intersect grid BBOX, minx=" + xCoords[0] + ", miny=" + yCoords[0]
+ ", maxx=" + xCoords[xCoords.length - 1] + ", maxy=" + yCoords[yCoords.length - 1]);
}
ProjectionRect.intersect(fullBB, projRect, projRect);
ProjectionPoint lowerLeft = projRect.getLowerLeftPoint();
ProjectionPoint upperRigth = projRect.getUpperRightPoint();
double minX = lowerLeft.getX();
double minY = lowerLeft.getY();
double maxX = upperRigth.getX();
double maxY = upperRigth.getY();
int minY_idx = yAxis.findCoordElement(minY);
int maxY_idx = yAxis.findCoordElement(maxY);
// When yAxis is increasing, "minY_idx < maxY_idx == true".
// When yAxis is decreasing, "minY_idx > maxY_idx == true".
// To handle both cases, we need to use the min() and max() functions.
Range yRange = new Range(Math.min(minY_idx, maxY_idx), Math.max(minY_idx, maxY_idx), horizStride);
int minX_idx = xAxis.findCoordElement(minX);
int maxX_idx = xAxis.findCoordElement(maxX);
Range xRange = new Range(Math.min(minX_idx, maxX_idx), Math.max(minX_idx, maxX_idx), horizStride);
yxRanges.add(0, yRange);
yxRanges.add(1, xRange);
}
private void addGlobalAttributes(NetcdfFileWriter writer, ucar.nc2.dt.GridDataset gds, LatLonRect llbb) {
// global attributes
for (Attribute att : gds.getGlobalAttributes()) {
if (att.getShortName().equals(CDM.FILE_FORMAT))
continue;
if (att.getShortName().equals(_Coordinate._CoordSysBuilder))
continue;
writer.addGroupAttribute(null, att);
}
Attribute att = gds.findGlobalAttributeIgnoreCase(CDM.CONVENTIONS);
if (att == null || !att.getStringValue().startsWith("CF-")) // preserve previous version of CF Convention if it
// exists
writer.addGroupAttribute(null, new Attribute(CDM.CONVENTIONS, "CF-1.0"));
writer.addGroupAttribute(null,
new Attribute("History", "Translated to CF-1.0 Conventions by Netcdf-Java CDM (CFGridWriter2)\n"
+ "Original Dataset = " + gds.getLocation() + "; Translation Date = " + CalendarDate.present()));
// this will replace any existing
writer.addGroupAttribute(null, new Attribute(ACDD.LAT_MIN, llbb.getLatMin()));
writer.addGroupAttribute(null, new Attribute(ACDD.LAT_MAX, llbb.getLatMax()));
writer.addGroupAttribute(null, new Attribute(ACDD.LON_MIN, llbb.getLonMin()));
writer.addGroupAttribute(null, new Attribute(ACDD.LON_MAX, llbb.getLonMax()));
}
private void addCFAnnotations(NetcdfFileWriter writer, ucar.nc2.dt.GridDataset gds, List gridList,
NetcdfDataset ncd, List axisList, boolean addLatLon) {
// Group root = ncfile.getRootGroup();
for (String gridName : gridList) {
GridDatatype grid = gds.findGridDatatype(gridName);
Variable newV = writer.findVariable(gridName);
if (newV == null) {
log.debug("NetcdfCFWriter cant find " + gridName + " in gds " + gds.getLocation());
continue;
}
// annotate Variable for CF
StringBuilder sbuff = new StringBuilder();
GridCoordSystem gcs = grid.getCoordinateSystem();
for (Variable axis : gcs.getCoordinateAxes()) {
sbuff.append(axis.getFullName()).append(" ");
}
if (addLatLon)
sbuff.append("lat lon");
newV.addAttribute(new Attribute(CF.COORDINATES, sbuff.toString()));
// looking for coordinate transform variables
for (CoordinateTransform ct : gcs.getCoordinateTransforms()) {
Variable v = ncd.findVariable(ct.getName());
if (ct.getTransformType() == TransformType.Projection)
newV.addAttribute(new Attribute(CF.GRID_MAPPING, v.getFullName()));
}
}
for (CoordinateAxis axis : axisList) {
Variable newV = writer.findVariable(axis.getFullNameEscaped());
if ((axis.getAxisType() == AxisType.Height) || (axis.getAxisType() == AxisType.Pressure)
|| (axis.getAxisType() == AxisType.GeoZ)) {
if (null != axis.getPositive())
newV.addAttribute(new Attribute(CF.POSITIVE, axis.getPositive()));
}
if (axis.getAxisType() == AxisType.Lat) {
newV.addAttribute(new Attribute(CDM.UNITS, "degrees_north"));
newV.addAttribute(new Attribute(CF.STANDARD_NAME, "latitude"));
}
if (axis.getAxisType() == AxisType.Lon) {
newV.addAttribute(new Attribute(CDM.UNITS, "degrees_east"));
newV.addAttribute(new Attribute(CF.STANDARD_NAME, "longitude"));
}
if (axis.getAxisType() == AxisType.GeoX) {
newV.addAttribute(new Attribute(CF.STANDARD_NAME, "projection_x_coordinate"));
}
if (axis.getAxisType() == AxisType.GeoY) {
newV.addAttribute(new Attribute(CF.STANDARD_NAME, "projection_y_coordinate"));
}
}
// coordinate transform variables : must convert false easting, northing to km
List ctvList = new ArrayList<>();
for (ucar.nc2.dt.GridDataset.Gridset gridSet : gds.getGridsets()) {
ucar.nc2.dt.GridCoordSystem gcs = gridSet.getGeoCoordSystem();
ProjectionCT pct = gcs.getProjectionCT();
if (pct != null) {
Variable v = writer.findVariable(pct.getName()); // look for the ctv
if ((v != null) && !ctvList.contains(v)) {
convertProjectionCTV((NetcdfDataset) gds.getNetcdfFile(), v.attributes());
ctvList.add(v);
}
}
}
}
private void convertProjectionCTV(NetcdfDataset ds, AttributeContainer ctv) {
Attribute att = ctv.findAttribute(_Coordinate.TransformType);
if ((null != att) && att.getStringValue().equals("Projection")) {
Attribute east = ctv.findAttribute("false_easting");
Attribute north = ctv.findAttribute("false_northing");
if ((null != east) || (null != north)) {
double scalef = AbstractTransformBuilder.getFalseEastingScaleFactor(ds, ctv);
if (scalef != 1.0) {
convertAttribute(ctv, east, scalef);
convertAttribute(ctv, north, scalef);
}
}
}
}
private void convertAttribute(AttributeContainer ctv, Attribute att, double scalef) {
if (att == null)
return;
double val = scalef * att.getNumericValue().doubleValue();
ctv.addAttribute(new Attribute(att.getShortName(), val));
}
private long addLatLon2D(NetcdfFile ncfile, List varList, Projection proj, CoordinateAxis xaxis,
CoordinateAxis yaxis) throws IOException {
double[] xData = (double[]) xaxis.read().get1DJavaArray(double.class);
double[] yData = (double[]) yaxis.read().get1DJavaArray(double.class);
List dims = new ArrayList<>();
dims.add(yaxis.getDimension(0));
dims.add(xaxis.getDimension(0));
Variable latVar = new Variable(ncfile, null, null, "lat");
latVar.setDataType(DataType.DOUBLE);
latVar.setDimensions(dims);
latVar.addAttribute(new Attribute(CDM.UNITS, CDM.LAT_UNITS));
latVar.addAttribute(new Attribute(CDM.LONG_NAME, "latitude coordinate"));
latVar.addAttribute(new Attribute(CF.STANDARD_NAME, "latitude"));
latVar.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.Lat.toString()));
Variable lonVar = new Variable(ncfile, null, null, "lon");
lonVar.setDataType(DataType.DOUBLE);
lonVar.setDimensions(dims);
lonVar.addAttribute(new Attribute(CDM.UNITS, CDM.LON_UNITS));
lonVar.addAttribute(new Attribute(CDM.LONG_NAME, "longitude coordinate"));
lonVar.addAttribute(new Attribute(CF.STANDARD_NAME, "longitude"));
lonVar.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.Lon.toString()));
int nx = xData.length;
int ny = yData.length;
// create the data
double[] latData = new double[nx * ny];
double[] lonData = new double[nx * ny];
for (int i = 0; i < ny; i++) {
for (int j = 0; j < nx; j++) {
ProjectionPoint projPoint = ProjectionPoint.create(xData[j], yData[i]);
LatLonPoint latlonPoint = proj.projToLatLon(projPoint);
latData[i * nx + j] = latlonPoint.getLatitude();
lonData[i * nx + j] = latlonPoint.getLongitude();
}
}
Array latDataArray = Array.factory(DataType.DOUBLE, new int[] {ny, nx}, latData);
latVar.setCachedData(latDataArray, false);
Array lonDataArray = Array.factory(DataType.DOUBLE, new int[] {ny, nx}, lonData);
lonVar.setCachedData(lonDataArray, false);
varList.add(latVar);
varList.add(lonVar);
long size = 0;
size += latVar.getSize() * latVar.getElementSize();
size += lonVar.getSize() * lonVar.getElementSize();
return size;
}
private Range makeVerticalRange(Range zRange, CoordinateAxis1D vertAxis) {
return (zRange != null) && (vertAxis != null) && (vertAxis.getSize() > 1) ? zRange : null;
}
private Range makeTimeRange(CalendarDateRange dateRange, CoordinateAxis1DTime timeAxis, int stride_time)
throws InvalidRangeException {
Range timeRange = null;
if ((dateRange != null) && (timeAxis != null)) {
int startIndex = timeAxis.findTimeIndexFromCalendarDate(dateRange.getStart());
int endIndex = timeAxis.findTimeIndexFromCalendarDate(dateRange.getEnd());
if (startIndex < 0)
throw new InvalidRangeException(
"start time=" + dateRange.getStart() + " must be >= " + timeAxis.getCalendarDate(0));
if (endIndex < 0)
throw new InvalidRangeException(
"end time=" + dateRange.getEnd() + " must be >= " + timeAxis.getCalendarDate(0));
if (stride_time <= 1)
stride_time = 1;
timeRange = new Range(startIndex, endIndex, stride_time);
}
return timeRange;
}
private void addCoordinateAxis(GridCoordSystem gcs, List varNameList, List varList,
List axisList) {
for (CoordinateAxis axis : gcs.getCoordinateAxes()) {
if (!varNameList.contains(axis.getFullName())) {
varNameList.add(axis.getFullName());
varList.add(axis); // LOOK axis hasnt been subset yet !!
axisList.add(axis);
// if (timeAxis != null && timeAxis.isInterval()) {
// LOOK gotta add the bounds !!!
// }
}
}
}
private void addCoordinateTransform(GridCoordSystem gcs, NetcdfFile ncd, List varNameList,
List varList) {
List coordTransformNames = new ArrayList<>();
// look for variables detected as coordinate transforms
for (CoordinateTransform ct : gcs.getCoordinateTransforms()) {
coordTransformNames.add(ct.getName());
}
// At the very least, LatLon_Projection will not show up in the coordinate transforms,
// so inspect the variables in varNameList for a grid_mapping attribute, and make sure it
// gets picked up as well
// see https://github.com/Unidata/python-workshop/issues/372
for (String varName : varNameList) {
Variable v = ncd.findVariable(varName);
String gridMappingName = v.attributes().findAttributeString(CF.GRID_MAPPING, null);
if (gridMappingName != null && (!coordTransformNames.contains(gridMappingName))) {
coordTransformNames.add(gridMappingName);
}
}
// add the coordinate transforms that we found
for (String ctVarName : coordTransformNames) {
Variable v = ncd.findVariable(ctVarName);
if (!varNameList.contains(ctVarName) && (null != v)) {
varNameList.add(ctVarName);
varList.add(v);
}
}
}
/*
* Process the coordinate transformations (formula_terms) and adds the variables needed for performing that
* transformation.
* Subsets the grids variables, if needed.
* Return size of variables added
*/
private long processTransformationVars(List varList, List varNameList, NetcdfDataset ncd,
ucar.nc2.dt.GridDataset gds, GridDatatype grid, Range timeRange, Range zRangeUse, Range yRange, Range xRange,
int y_stride, int x_stride) throws InvalidRangeException {
long varsSize = 0L;
List cctt = grid.getCoordinateSystem().getCoordinateTransforms();
for (CoordinateTransform ct : cctt) {
Parameter param = ct.findParameterIgnoreCase(CF.FORMULA_TERMS);
if (param != null) {
String[] varStrings = param.getStringValue().split(" ");
for (int i = 1; i < varStrings.length; i += 2) {
Variable paramVar = ncd.findVariable(varStrings[i].trim());
if (!varNameList.contains(varStrings[i]) && (null != paramVar)) {
if (gds.findGridDatatype(paramVar.getFullName()) != null) {
// Subset if needed
if ((null != timeRange) || (zRangeUse != null) || (x_stride > 1 && y_stride > 1)
|| (yRange != null || xRange != null)) {
GridDatatype complementaryGrid = gds.findGridDatatype(paramVar.getFullName());
complementaryGrid = complementaryGrid.makeSubset(null, null, timeRange, zRangeUse, yRange, xRange);
paramVar = complementaryGrid.getVariable();
}
} else {
// Also have to subset the var if it is not a grid but has vertical dimension (the dimensionless vars in
// the formula) and zRangeUse != null
if (zRangeUse != null && paramVar.getRank() == 1) {
List ranges = new ArrayList<>();
ranges.add(zRangeUse);
paramVar = paramVar.section(ranges);
}
}
varNameList.add(paramVar.getFullName());
varsSize += paramVar.getSize() * paramVar.getElementSize();
varList.add(paramVar);
}
}
}
}
return varsSize;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy