ucar.nc2.geotiff.GeoTiffWriter2 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!
package ucar.nc2.geotiff;
import ucar.ma2.*;
import ucar.nc2.Attribute;
import ucar.nc2.dataset.CoordinateAxis1D;
import ucar.nc2.dataset.CoordinateAxis2D;
import ucar.nc2.dt.GridCoordSystem;
import ucar.nc2.dt.GridDataset;
import ucar.nc2.dt.GridDatatype;
import ucar.unidata.geoloc.*;
import java.io.IOException;
/**
* Not sure where this is used - IDV ??
*
* @author caron
* @since 3/15/13
*/
public class GeoTiffWriter2 extends GeotiffWriter {
public GeoTiffWriter2(String fileOut) {
super(fileOut);
}
public void writeGrid(String gridDataset_filename, String gridName, int time, int level, boolean greyScale, LatLonRect pt) throws IOException {
double scaler;
try (GridDataset dataset = ucar.nc2.dt.grid.GridDataset.open(gridDataset_filename)) {
GridDatatype grid = dataset.findGridDatatype(gridName);
if (grid == null) {
throw new IllegalArgumentException("No grid named " + gridName + " in fileName");
}
GridCoordSystem gcs = grid.getCoordinateSystem();
ProjectionImpl proj = grid.getProjection();
if (!gcs.isRegularSpatial()) {
Attribute att = dataset.findGlobalAttributeIgnoreCase("datasetId");
if (att != null && att.getStringValue().contains("DMSP")) { // LOOK!!
writeSwathGrid(gridDataset_filename, gridName, time, level, greyScale, pt);
return;
} else {
throw new IllegalArgumentException("Must have 1D x and y axes for " + grid.getFullName());
}
}
CoordinateAxis1D xaxis = (CoordinateAxis1D) gcs.getXHorizAxis();
CoordinateAxis1D yaxis = (CoordinateAxis1D) gcs.getYHorizAxis();
if (!xaxis.isRegular() || !yaxis.isRegular()) {
throw new IllegalArgumentException("Must be evenly spaced grid = " + grid.getFullName());
}
// read in data
Array data = grid.readDataSlice(time, level, -1, -1);
Array lon = xaxis.read();
Array lat = yaxis.read();
//latlon coord does not need to time 1000.0
if (gcs.isLatLon()) {
scaler = 1.0;
} else {
scaler = 1000.0;
}
if (yaxis.getCoordValue(0) < yaxis.getCoordValue(1)) {
data = data.flip(0);
lat = lat.flip(0);
}
if (gcs.isLatLon()) {
data = geoShiftDataAtLon(data, lon);
lon = geoShiftLon(lon);
}
// now it is time to subset the data out of latlonrect
// it is assumed that latlonrect pt is in +-180
LatLonPointImpl llp0 = pt.getLowerLeftPoint();
LatLonPointImpl llpn = pt.getUpperRightPoint();
double minLon = llp0.getLongitude();
double minLat = llp0.getLatitude();
double maxLon = llpn.getLongitude();
double maxLat = llpn.getLatitude();
// (x1, y1) is upper left point and (x2, y2) is lower right point
int x1;
int x2;
int y1;
int y2;
double xStart;
double yStart;
if (!gcs.isLatLon()) {
ProjectionPoint pjp0 = proj.latLonToProj(maxLat, minLon);
x1 = getXIndex(lon, pjp0.getX(), 0);
y1 = getYIndex(lat, pjp0.getY(), 0);
yStart = pjp0.getY() * 1000.0; //latArray[y1];
xStart = pjp0.getX() * 1000.0; //lonArray[x1];
ProjectionPoint pjpn = proj.latLonToProj(minLat, maxLon);
x2 = getXIndex(lon, pjpn.getX(), 1);
y2 = getYIndex(lat, pjpn.getY(), 1);
} else {
xStart = minLon;
yStart = maxLat;
x1 = getLonIndex(lon, minLon, 0);
y1 = getLatIndex(lat, maxLat, 0);
x2 = getLonIndex(lon, maxLon, 1);
y2 = getLatIndex(lat, minLat, 1);
}
// data must go from top to bottom
double xInc = xaxis.getIncrement() * scaler;
double yInc = Math.abs(yaxis.getIncrement()) * scaler;
// subsetting data inside the box
Array data1 = getYXDataInBox(data, x1, x2, y1, y2);
if (pageNumber > 1) {
geotiff.initTags();
}
// write it out
writeGrid(grid, data1, greyScale, xStart, yStart, xInc, yInc, pageNumber);
pageNumber++;
}
}
/**
* Write Swath Grid data to the geotiff file.
*
* @param fileName _more_
* @param gridName _more_
* @param time _more_
* @param level _more_
* @param greyScale _more_
* @param llr _more_
* @throws IOException _more_
*/
void writeSwathGrid(String fileName, String gridName, int time, int level, boolean greyScale, LatLonRect llr) throws IOException {
double scaler;
GridDataset dataset = ucar.nc2.dt.grid.GridDataset.open(fileName);
GridDatatype grid = dataset.findGridDatatype(gridName);
GridCoordSystem gcs = grid.getCoordinateSystem();
ProjectionImpl proj = grid.getProjection();
CoordinateAxis2D xaxis = (CoordinateAxis2D) gcs.getXHorizAxis();
CoordinateAxis2D yaxis = (CoordinateAxis2D) gcs.getYHorizAxis();
// read in data
Array data = grid.readDataSlice(time, level, -1, -1);
Array lon = xaxis.read();
Array lat = yaxis.read();
double[] swathInfo = getSwathLatLonInformation(lat, lon);
//latlon coord does not need to time 1000.0
if (gcs.isLatLon()) {
scaler = 1.0;
} else {
scaler = 1000.0;
}
//if (yaxis.getCoordValue(0, 0) < yaxis.getCoordValue(0, 1)) {//???
data = data.flip(0);
//lat = lat.flip(0);
//}
if (gcs.isLatLon()) {
data = geoShiftDataAtLon(data, lon);
lon = geoShiftLon(lon);
}
double minLon;
double minLat;
double maxLon;
double maxLat;
double xStart;
double yStart; //upper right point
double xInc = swathInfo[0] * scaler;
double yInc = swathInfo[1] * scaler;
// (x1, y1) is upper left point and (x2, y2) is lower right point
int x1;
int x2;
int y1;
int y2;
if (llr == null) { //get the whole area
minLon = swathInfo[4];
minLat = swathInfo[2];
maxLon = swathInfo[5];
maxLat = swathInfo[3];
xStart = minLon;
yStart = maxLat;
x1 = 0;
y1 = 0;
x2 = (int) ((maxLon - minLon) / xInc + 0.5);
y2 = (int) ((maxLat - minLat) / yInc + 0.5);
} else { //assign the special area surrounded by the llr
LatLonPointImpl llp0 = llr.getLowerLeftPoint();
LatLonPointImpl llpn = llr.getUpperRightPoint();
minLon = (llp0.getLongitude() < swathInfo[4])
? swathInfo[4]
: llp0.getLongitude();
minLat = (llp0.getLatitude() < swathInfo[2])
? swathInfo[2]
: llp0.getLatitude();
maxLon = (llpn.getLongitude() > swathInfo[5])
? swathInfo[5]
: llpn.getLongitude();
maxLat = (llpn.getLatitude() > swathInfo[3])
? swathInfo[3]
: llpn.getLatitude();
//construct the swath LatLonRect
LatLonPointImpl pUpLeft = new LatLonPointImpl(swathInfo[3], swathInfo[4]);
LatLonPointImpl pDownRight = new LatLonPointImpl(swathInfo[2], swathInfo[5]);
LatLonRect swathLLR = new LatLonRect(pUpLeft, pDownRight);
LatLonRect bIntersect = swathLLR.intersect(llr);
if (bIntersect == null) {
throw new IllegalArgumentException("The assigned extent of latitude and longitude is unvalid. No intersection with the swath extent");
}
xStart = minLon;
yStart = maxLat;
x1 = (int) ((minLon - swathInfo[4]) / xInc + 0.5);
y1 = (int) Math.abs((maxLat - swathInfo[3]) / yInc + 0.5);
x2 = (int) ((maxLon - swathInfo[4]) / xInc + 0.5);
y2 = (int) Math.abs((minLat - swathInfo[3]) / yInc + 0.5);
}
if (!gcs.isLatLon()) {
ProjectionPoint pjp0 = proj.latLonToProj(maxLat, minLon);
x1 = getXIndex(lon, pjp0.getX(), 0);
y1 = getYIndex(lat, pjp0.getY(), 0);
yStart = pjp0.getY() * 1000.0; //latArray[y1];
xStart = pjp0.getX() * 1000.0; //lonArray[x1];
ProjectionPoint pjpn = proj.latLonToProj(minLat, maxLon);
x2 = getXIndex(lon, pjpn.getX(), 1);
y2 = getYIndex(lat, pjpn.getY(), 1);
} else {
//calculate the x1, x2, y1, y2, xstart, ystart.
}
Array targetImage = getTargetImagerFromSwath(lat, lon, data, swathInfo);
Array interpolatedImage = interpolation(targetImage);
Array clippedImage = getClippedImageFromInterpolation(interpolatedImage, x1, x2, y1, y2);
//Array clippedImage = getYXDataInBox(interpolatedImage, x1, x2, y1, y2);
if (pageNumber > 1) {
geotiff.initTags();
}
writeGrid(grid, clippedImage, greyScale, xStart, yStart, xInc, yInc, pageNumber);
pageNumber++;
}
int getXIndex(Array aAxis, double value, int side) {
IndexIterator aIter = aAxis.getIndexIterator();
int count = 0;
int isInd = 0;
double aValue = aIter.getFloatNext();
if ((aValue == value) || (aValue > value)) {
return 0;
}
while (aIter.hasNext() && (aValue < value)) {
count++;
aValue = aIter.getFloatNext();
if (aValue == value) {
isInd = 1;
}
}
if (isInd == 1) {
count += side;
}
count -= side;
return count;
}
int getYIndex(Array aAxis, double value, int side) {
IndexIterator aIter = aAxis.getIndexIterator();
int count = 0;
int isInd = 0;
double aValue = aIter.getFloatNext();
if ((aValue == value) || (aValue < value)) {
return 0;
}
while (aIter.hasNext() && (aValue > value)) {
count++;
aValue = aIter.getFloatNext();
if (aValue == value) {
isInd = 1;
}
}
if (isInd == 1) {
count += side;
}
count -= side;
return count;
}
int getLatIndex(Array lat, double value, int side) {
IndexIterator latIter = lat.getIndexIterator();
int count = 0;
int isInd = 0;
//LatLonPoint p0 = new LatLonPointImpl(lat.getFloat(ind.set(0)), 0);
double xlat = latIter.getFloatNext();
if (xlat == value) {
return 0;
}
while (latIter.hasNext() && (xlat > value)) {
count++;
xlat = latIter.getFloatNext();
if (xlat == value) {
isInd = 1;
}
}
if (isInd == 1) {
count += side;
}
count -= side;
return count;
}
int getLonIndex(Array lon, double value, int side) {
IndexIterator lonIter = lon.getIndexIterator();
int count = 0;
int isInd = 0;
// double xlon = lon.getFloat(ind.set(0));
float xlon = lonIter.getFloatNext();
if (xlon > 180) {
xlon = xlon - 360;
}
if (xlon == value) {
return 0;
}
while (lonIter.hasNext() && (xlon < value)) {
count++;
xlon = lonIter.getFloatNext();
if (xlon > 180) {
xlon = xlon - 360;
}
if (xlon == value) {
isInd = 1;
}
}
if (isInd == 1) {
count += side;
}
count -= side;
return count;
}
Array getYXDataInBox(Array data, int x1, int x2, int y1,
int y2) throws java.io.IOException {
int rank = data.getRank();
int[] start = new int[rank];
int[] shape = new int[rank];
for (int i = 0; i < rank; i++) {
start[i] = 0;
shape[i] = 1;
}
if ((y1 >= 0) && (y2 >= 0)) {
start[0] = y1;
shape[0] = y2 - y1;
}
if ((x1 >= 0) && (x2 >= 0)) {
start[1] = x1;
shape[1] = x2 - x1;
}
// read it
Array dataVolume;
try {
dataVolume = data.section(start, shape);
} catch (Exception e) {
throw new java.io.IOException(e.getMessage());
}
return dataVolume;
}
Array getClippedImageFromInterpolation(Array arr, int x1, int x2, int y1, int y2) {
int[] srcShape = arr.getShape();
int rank = arr.getRank();
int[] start = new int[rank];
int[] shape = new int[rank];
for (int i = 0; i < rank; i++) {
start[i] = 0;
shape[i] = 1;
}
if ((y1 >= 0) && (y2 >= 0)) {
start[0] = y1;
shape[0] = y2 - y1;
}
if ((x1 >= 0) && (x2 >= 0)) {
start[1] = x1;
shape[1] = x2 - x1;
}
Array dataVolume = Array.factory(DataType.FLOAT, shape);
int count = 0;
for (int i = y1; i < y2; i++) {
for (int j = x1; j < x2; j++) {
int index = i * srcShape[1] + j;
if (index >= srcShape[0] * srcShape[1]) {
index = srcShape[0] * srcShape[1] - 1;
}
float curValue = arr.getFloat(index);
if (count >= shape[0] * shape[1]) {
count = shape[0] * shape[1] - 1;
}
dataVolume.setFloat(count, curValue);
count++;
}
}
return dataVolume;
}
Array getTargetImagerFromSwath(Array lat, Array lon, Array data, double[] swathInfo) {
int srcDataHeight = data.getShape()[0];
int srcDataWidth = data.getShape()[1];
int BBoxHeight = (int) ((swathInfo[3] - swathInfo[2]) / swathInfo[1] + 0.5);
int BBoxWidth = (int) ((swathInfo[5] - swathInfo[4]) / swathInfo[0] + 0.5);
int BBoxShape[] = {BBoxHeight, BBoxWidth}; //[height, width]
Array bBoxArray = Array.factory(DataType.FLOAT, BBoxShape);
double startLon, startLat; //upper left and lower right
startLon = swathInfo[4];
startLat = swathInfo[3];
IndexIterator dataIter = data.getIndexIterator();
IndexIterator latIter = lat.getIndexIterator();
IndexIterator lonIter = lon.getIndexIterator();
for (int i = 0; i < srcDataHeight; i++) {
for (int j = 0; j < srcDataWidth; j++) {
while (latIter.hasNext() && lonIter.hasNext()
&& dataIter.hasNext()) {
float curLat = latIter.getFloatNext();
float curLon = lonIter.getFloatNext();
float curPix = dataIter.getFloatNext();
float alreadyValue;
int curPixelInBBoxIndex =
getIndexOfBBFromLatlonOfOri(startLat, startLon,
swathInfo[1], swathInfo[0], curLat, curLon,
BBoxHeight, BBoxWidth);
try {
alreadyValue = bBoxArray.getFloat(curPixelInBBoxIndex);
} catch (Exception e) {
alreadyValue = 0;
}
if (alreadyValue > 0) { //This pixel had been filled. So calculate the average
bBoxArray.setFloat(curPixelInBBoxIndex, (curPix + alreadyValue) / 2);
} else {
bBoxArray.setFloat(curPixelInBBoxIndex, curPix);
}
}
}
}
return bBoxArray;
}
int getIndexOfBBFromLatlonOfOri(double sLat, double sLon, //LatLon of the start point
double latInc, double lonInc, //The increment in Lat/Lon direction
double curLat, double curLon, //The current Lat/Lon read from the swath image
int bbHeight, int bbWidth) //The width and height of target image
{
double lonDelta = Math.abs((curLon - sLon) / lonInc);
double latDelta = Math.abs((curLat - sLat) / latInc);
int row = (int) (lonDelta + 0.5);
if (row >= bbWidth - 1) {
row = bbWidth - 2;
}
if (row == 0) {
row = 1;
}
int col = (int) (latDelta + 0.5);
int index = col * bbWidth + row;
if (index >= bbHeight * bbWidth) {
index = (col - 1) * bbWidth + row;
}
return index;
}
/**
* interpolate the swath data to regular grid
*
* @param arr swath array
* @return regular grid
*/
Array interpolation(Array arr) {
int[] orishape = arr.getShape();
int width = orishape[1];
int height = orishape[0];
int pixelNum = width * height;
Array interpolatedArray = Array.factory(DataType.FLOAT, orishape);
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
int curIndex = i * width + j;
float curValue = arr.getFloat(curIndex);
if (curValue == 0) //Black hole. Need to fill.
{
float tempPixelSum = 0;
int numNeighborHasValue = 0;
//Get the values of eight neighborhood
if ((curIndex - 1 >= 0) && (curIndex - 1 < pixelNum)) {
float left = arr.getFloat(curIndex - 1);
if (left > 0) {
tempPixelSum += left;
numNeighborHasValue++;
}
}
if ((curIndex + 1 >= 0) && (curIndex + 1 < pixelNum)) {
float right = arr.getFloat(curIndex + 1);
if (right > 0) {
tempPixelSum += right;
numNeighborHasValue++;
}
}
if ((curIndex - width >= 0) && (curIndex - width < pixelNum)) {
float up = arr.getFloat(curIndex - width);
if (up > 0) {
tempPixelSum += up;
numNeighborHasValue++;
}
}
if ((curIndex + width >= 0) && (curIndex + width < pixelNum)) {
float down = arr.getFloat(curIndex + width);
if (down > 0) {
tempPixelSum += down;
numNeighborHasValue++;
}
}
if ((curIndex - width - 1 >= 0) && (curIndex - width - 1 < pixelNum)) {
float upleft = arr.getFloat(curIndex - width - 1);
if (upleft > 0) {
tempPixelSum += upleft;
numNeighborHasValue++;
}
}
if ((curIndex - width + 1 >= 0) && (curIndex - width + 1 < pixelNum)) {
float upright = arr.getFloat(curIndex - width + 1);
if (upright > 0) {
tempPixelSum += upright;
numNeighborHasValue++;
}
}
if ((curIndex + width - 1 >= 0) && (curIndex + width - 1 < pixelNum)) {
float downleft = arr.getFloat(curIndex + width - 1);
if (downleft > 0) {
tempPixelSum += downleft;
numNeighborHasValue++;
}
}
if ((curIndex + width + 1 >= 0) && (curIndex + width + 1 < pixelNum)) {
float downright = arr.getFloat(curIndex + width + 1);
if (downright > 0) {
tempPixelSum += downright;
numNeighborHasValue++;
}
}
if (tempPixelSum > 0) {
float val = numNeighborHasValue == 0 ? 0 : tempPixelSum / numNeighborHasValue;
interpolatedArray.setFloat(curIndex, val);
}
} else {
interpolatedArray.setFloat(curIndex, curValue);
}
}
}
return interpolatedArray;
}
/**
* get lat lon information from the swath
*
* @param lat 2D lat array
* @param lon 2D lon array
* @return double array for geotiff
*/
double[] getSwathLatLonInformation(Array lat, Array lon) {
// Calculate the increment of latitude and longitude of original swath data
// Calculate the size of the boundingBox
// element0: Longitude increment
// element1: Latitude increment
// element2: minLat
// element3: maxLat
// element4: minLon
// element5: maxLon
// element6: width of the boundingBox
// element7: height of the boundingBox
double increment[] = {0, 0, 0, 0, 0, 0, 0, 0};
IndexIterator latIter = lat.getIndexIterator();
IndexIterator lonIter = lon.getIndexIterator();
int numScan = (lat.getShape())[0];
int numSample = (lat.getShape())[1];
double maxLat = -91,
minLat = 91,
maxLon = -181,
minLon = 181;
float firstLineStartLat = 0;
float firstLineStartLon = 0;
float firstLineEndLat = 0;
float firstLineEndLon = 0;
float lastLineStartLat = 0;
float lastLineStartLon = 0;
float lastLineEndLat = 0;
float lastLineEndLon = 0;
for (int i = 0; i < numScan; i++) {
for (int j = 0; j < numSample; j++) {
if (latIter.hasNext() && lonIter.hasNext()) {
float curLat = latIter.getFloatNext();
float curLon = lonIter.getFloatNext();
if ((i == 0) && (j == 0)) {
firstLineStartLat = curLat;
firstLineStartLon = curLon;
} else if ((i == 0) && (j == numSample - 1)) {
firstLineEndLat = curLat;
firstLineEndLon = curLon;
} else if ((i == numScan - 1) && (j == 0)) {
lastLineStartLat = curLat;
lastLineStartLon = curLon;
} else if ((i == numScan - 1) && (j == numSample - 1)) {
lastLineEndLat = curLat;
lastLineEndLon = curLon;
}
}
}
}
double[] edgeLat = {firstLineStartLat, firstLineEndLat, lastLineStartLat, lastLineEndLat};
double[] edgeLon = {firstLineStartLon, firstLineEndLon, lastLineStartLon, lastLineEndLon};
for (int i = 0; i < edgeLat.length; i++) {
maxLat = ((maxLat > edgeLat[i])
? maxLat
: edgeLat[i]);
minLat = ((minLat < edgeLat[i])
? minLat
: edgeLat[i]);
maxLon = ((maxLon > edgeLon[i])
? maxLon
: edgeLon[i]);
minLon = ((minLon < edgeLon[i])
? minLon
: edgeLon[i]);
}
double xInc1 = Math.abs((firstLineEndLon - firstLineStartLon) / numSample);
//double xInc2 = Math.abs((lastLineEndLon - lastLineStartLon)/numSample);
double yInc1 = Math.abs((lastLineStartLat - firstLineStartLat) / numScan);
//double yInc2 = Math.abs((lastLineEndLat - firstLineEndLat)/numScan);
increment[0] = xInc1; // > xInc2 ? xInc1 : xInc2;
increment[1] = yInc1; // > yInc2 ? yInc1 : yInc2;
increment[2] = minLat;
increment[3] = maxLat;
increment[4] = minLon;
increment[5] = maxLon;
increment[6] = (maxLon - minLon) / xInc1;
increment[7] = (maxLat - minLat) / yInc1;
return increment;
}
private Array geoShiftLon(Array lon) {
int count = 0;
Index lonIndex = lon.getIndex();
int[] lonShape = lon.getShape();
ArrayFloat slon = new ArrayFloat(new int[]{lonShape[0]});
Index slonIndex = slon.getIndex();
IndexIterator lonIter = lon.getIndexIterator();
LatLonPoint p0 = new LatLonPointImpl(0, lon.getFloat(lonIndex.set(lonShape[0] - 1)));
LatLonPoint pN = new LatLonPointImpl(0, lon.getFloat(lonIndex.set(0)));
while (lonIter.hasNext()) {
float l = lonIter.getFloatNext();
if (l > 180.0) {
count++;
}
}
//checking if the 0 point and the N point are the same point
int spoint;
if (p0.getLongitude() == pN.getLongitude()) {
spoint = lonShape[0] - count - 1;
} else {
spoint = lonShape[0] - count;
}
if ((count > 0) && (lonShape[0] > count)) {
for (int j = 1; j < lonShape[0]; j++) {
int jj;
if (j >= count) {
jj = j - count;
} else {
jj = j + spoint;
}
float dd = lon.getFloat(lonIndex.set(jj));
slon.setFloat(slonIndex.set(j), (float) LatLonPointImpl.lonNormal(dd));
}
if (p0.getLongitude() == pN.getLongitude()) {
float dd = slon.getFloat(slonIndex.set(lonShape[0] - 1));
slon.setFloat(slonIndex.set(0), -(float) LatLonPointImpl.lonNormal(dd));
}
return slon;
} else {
return lon;
}
}
}
© 2015 - 2024 Weber Informatics LLC | Privacy Policy