Issue #2190 Sped up viirs map projection for use with nucaps cloud height

Change-Id: I7a99e734b9aede618f3afe7e01f3103f92c3fa52

Former-commit-id: 5169afc30f [formerly 3bc32de5b1029287403173cc97d815e018a4531e]
Former-commit-id: 8c26afa091
This commit is contained in:
Max Schenkelberg 2013-08-23 16:23:04 -05:00
parent 01edf7207c
commit 3425c21404
19 changed files with 678 additions and 656 deletions

View file

@ -215,10 +215,15 @@ public class CloudHeightCalculatorPorted {
// solve for T of Parcel Trajectory at P
presTrajL = ((int) ((sounding.get(j).getPressure() + 49.0) / 50)) * 50;
presTrajU = presTrajL - 50;
Tparcel = (parcelTrajectory[20 - (int) (presTrajU / 50)] - parcelTrajectory[20 - (int) (presTrajL / 50)])
int lIdx = 20 - (int) (presTrajL / 50);
int uIdx = 20 - (int) (presTrajU / 50);
if (lIdx < 0 || uIdx >= parcelTrajectory.length) {
continue;
}
Tparcel = (parcelTrajectory[uIdx] - parcelTrajectory[lIdx])
/ -50
* (sounding.get(j).getPressure() - presTrajL)
+ parcelTrajectory[20 - (int) (presTrajL / 50)];
+ parcelTrajectory[lIdx];
if (sounding.get(j).getTemperature() > Tparcel) {
presEL = sounding.get(j - 1).getPressure();
break;

View file

@ -452,11 +452,6 @@ public class CloudHeightResource extends
int x = (int) c.x;
int y = (int) c.y;
int nx = spatialObject.getNx();
int ny = spatialObject.getNy();
if (x < 0 || x >= nx || y < 0 || y >= ny) {
return temps;
}
double temp = getTemperature(rsc, latLon);
if (Double.isNaN(temp)) {
@ -474,15 +469,10 @@ public class CloudHeightResource extends
for (j = yStart; j < yEnd; j++) {
jj = y + j;
if (jj < 0 || jj >= ny) {
continue;
}
for (i = xStart; i < xEnd; i++) {
ii = x + i;
if (ii >= 0 && ii < nx) {
elements[numElements++] = getTemperature(rsc,
spatialObject, ii, jj);
}
elements[numElements++] = getTemperature(rsc, spatialObject,
ii, jj);
}
}

View file

@ -23,7 +23,6 @@ import org.geotools.coverage.grid.GeneralGridGeometry;
import org.geotools.coverage.grid.GridEnvelope2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.referencing.operation.DefaultMathTransformFactory;
import org.opengis.coverage.grid.GridEnvelope;
import org.opengis.geometry.Envelope;
@ -34,7 +33,6 @@ import org.opengis.referencing.operation.TransformException;
import com.raytheon.uf.common.geospatial.CRSCache;
import com.raytheon.uf.common.geospatial.util.EnvelopeIntersection;
import com.raytheon.uf.common.geospatial.util.WorldWrapCorrector;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
@ -50,6 +48,7 @@ import com.vividsolutions.jts.geom.Geometry;
* ------------ ---------- ----------- --------------------------
* Aug 8, 2012 mschenke Initial creation
* Apr 03, 2013 1824 bsteffen Fix Tile Geometry creation.
* Aug 27, 2013 2190 mschenke Remove unused transforms
*
* </pre>
*
@ -72,12 +71,6 @@ public class TileLevel {
private MathTransform tileCRSToTargetGrid;
/** Target grid to lat/lon for world wrap correcting */
private MathTransform targetGridToLatLon;
/** World wrap corrector, corrects Tile borders */
private WorldWrapCorrector corrector;
/** Level of this TileLevel */
private int tileLevel;
@ -130,12 +123,6 @@ public class TileLevel {
targetGeometry.getCoordinateReferenceSystem()),
targetGeometry.getGridToCRS(PixelInCell.CELL_CORNER)
.inverse());
targetGridToLatLon = factory.createConcatenatedTransform(
targetGeometry.getGridToCRS(PixelInCell.CELL_CORNER),
CRSCache.getInstance().findMathTransform(
targetGeometry.getCoordinateReferenceSystem(),
DefaultGeographicCRS.WGS84));
corrector = new WorldWrapCorrector(targetGeometry);
Envelope levelEnv = levelGeometry.getEnvelope();
double[] in = new double[] {
@ -158,14 +145,18 @@ public class TileLevel {
+ targetGeometry.getGridRange().getSpan(1) * 0.75;
}
double[] input = new double[] { mapPointX, mapPointY,
double[] targetGridPoints = new double[] { mapPointX, mapPointY,
mapPointX + 1, mapPointY + 1 };
double[] output = new double[input.length];
double[] tileCRSPoints = new double[targetGridPoints.length];
double[] tileGridPoints = new double[tileCRSPoints.length];
tileCRSToTargetGrid.inverse().transform(input, 0, output, 0, 2);
crsToGrid.transform(output, 0, input, 0, 2);
pixelDensity = 1.0 / Math.abs(new Coordinate(input[0], input[1],
0.0).distance(new Coordinate(input[2], input[3], 0.0)));
tileCRSToTargetGrid.inverse().transform(targetGridPoints, 0,
tileCRSPoints, 0, 2);
crsToGrid.transform(tileCRSPoints, 0, tileGridPoints, 0, 2);
pixelDensity = 1.0 / Math.abs(new Coordinate(tileGridPoints[0],
tileGridPoints[1], 0.0).distance(new Coordinate(
tileGridPoints[2], tileGridPoints[3], 0.0)));
} catch (Exception e) {
throw new RuntimeException(
"Cannot tranform tile CRS into target CRS", e);
@ -285,9 +276,7 @@ public class TileLevel {
int tileWidth = Math.min(endX - tileX, tileSize);
int tileHeight = Math.min(endY - tileY, tileSize);
range = new GridEnvelope2D(tileX, tileY, tileWidth,
tileHeight);
range = new GridEnvelope2D(tileX, tileY, tileWidth, tileHeight);
Envelope envelope = null;
// Convert grid range into crs envelope range
try {

View file

@ -229,6 +229,25 @@ public class TileSetRenderable implements IRenderable {
}
}
/**
* Returns the {@link GeneralGridGeometry} the {@link TileSet} is currently
* projected for
*
* @return
*/
public GeneralGridGeometry getTargetGeometry() {
return tileSet != null ? tileSet.getTargetGeometry() : null;
}
/**
* Returns the {@link GridGeometry2D} of the {@link TileSet}
*
* @return
*/
public GridGeometry2D getTileSetGeometry() {
return tileSetGeometry;
}
/*
* (non-Javadoc)
*

View file

@ -48,6 +48,7 @@ import com.raytheon.uf.common.wxmath.ZToPsa;
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Aug 14, 2013 2260 bsteffen Initial creation
* Aug 27, 2013 2190 mschenke Fixed unit for VerticalSounding creation
*
* </pre>
*
@ -61,9 +62,9 @@ public class SoundingLayerBuilder {
private static final Unit<Length> HEIGHT_UNIT = SI.METER;
private static final Unit<Temperature> TEMPERATURE_UNIT = SI.CELSIUS;
private static final Unit<Temperature> TEMPERATURE_UNIT = SI.KELVIN;
private static final Unit<Temperature> DEWPOINT_UNIT = SI.CELSIUS;
private static final Unit<Temperature> DEWPOINT_UNIT = SI.KELVIN;
private static final Unit<Angle> WIND_DIRECTION_UNIT = NonSI.DEGREE_ANGLE;

View file

@ -64,6 +64,7 @@ import com.vividsolutions.jts.index.strtree.STRtree;
* ------------ ---------- ----------- --------------------------
* Jul 25, 2013 2190 mschenke Initial creation
* Aug 15, 2013 2260 bsteffen Switch poessounding to NSharp.
* Aug 27, 2013 2190 mschenke Fixed point query request
*
* </pre>
*
@ -77,14 +78,14 @@ public class NucapsSoundingProvider extends
private static final IUFStatusHandler statusHandler = UFStatus
.getHandler(NucapsSoundingProvider.class);
private static String[] SOUNDING_PARAMS = {
NucapsRecord.PDV_SURFACE_PRESSURE, NucapsRecord.PDV_PRESSURE,
NucapsRecord.PDV_TEMPERATURE,
private static String[] SOUNDING_PARAMS = { NucapsRecord.LATITUDE,
NucapsRecord.LONGITUDE, NucapsRecord.PDV_SURFACE_PRESSURE,
NucapsRecord.PDV_PRESSURE, NucapsRecord.PDV_TEMPERATURE,
NucapsRecord.PDV_WATER_VAPOR_MIXING_RATIO };
private static final double ENVELOPE_DISTANCE_DEG = 1.0;
private static final double MAX_SOUNDING_DISTANCE_METERS = 50 * 1000;
private static final double MAX_SOUNDING_DISTANCE_METERS = 100 * 1000;
private static final long GROUP_TIMERANGE_MILLIS = 15 * TimeUtil.MILLIS_PER_MINUTE;
@ -206,9 +207,10 @@ public class NucapsSoundingProvider extends
STRtree tree = new STRtree();
try {
PointDataContainer pdc = PointDataRequest.requestPointData(
time.getValidPeriod(), NucapsRecord.PLUGIN_NAME,
SOUNDING_PARAMS, null, constraints);
PointDataContainer pdc = PointDataRequest
.requestPointDataAllLevels(time.getValidPeriod(),
NucapsRecord.PLUGIN_NAME, SOUNDING_PARAMS, null,
constraints);
int size = pdc.getCurrentSz();
for (int i = 0; i < size; ++i) {
NucapsVerticalSounding sounding = new NucapsVerticalSounding(

View file

@ -2,4 +2,5 @@ source.. = src/
output.. = bin/
bin.includes = META-INF/,\
.,\
localization/
localization/,\
plugin.xml

View file

@ -0,0 +1,31 @@
<?xml version="1.0" encoding="UTF-8"?>
<!--
This_software_was_developed_and_/_or_modified_by_Raytheon_Company,
pursuant_to_Contract_DG133W-05-CQ-1067_with_the_US_Government.
U.S._EXPORT_CONTROLLED_TECHNICAL_DATA
This_software_product_contains_export-restricted_data_whose
export/transfer/disclosure_is_restricted_by_U.S._law._Dissemination
to_non-U.S._persons_whether_in_the_United_States_or_abroad_requires
an_export_license_or_other_authorization.
Contractor_Name:________Raytheon_Company
Contractor_Address:_____6825_Pine_Street,_Suite_340
________________________Mail_Stop_B8
________________________Omaha,_NE_68106
________________________402.291.0100
See_the_AWIPS_II_Master_Rights_File_("Master_Rights_File.pdf")_for
further_licensing_information.
-->
<?eclipse version="3.2"?>
<plugin>
<extension
point="com.raytheon.uf.viz.core.resource">
<resource
class="com.raytheon.uf.viz.npp.sounding.rsc.NPPSoundingMapResource"
name="NPP Sounding Availability"
renderingOrderId="PLOT"
resourceType="PLAN_VIEW"/>
</extension>
</plugin>

View file

@ -26,6 +26,7 @@ import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
@ -68,6 +69,7 @@ import com.vividsolutions.jts.geom.Coordinate;
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Jan 14, 2013 mschenke Initial creation
* Aug 27, 2013 2190 mschenke Fixed so can't click point when not visible
*
* </pre>
*
@ -80,7 +82,7 @@ public class NPPSoundingMapResource extends
private NPPSoundingMapInputHandler inputManager;
private Collection<NPPSoundingRecord> allRecords = new ArrayList<NPPSoundingRecord>();
private Collection<NPPSoundingRecord> allRecords = new HashSet<NPPSoundingRecord>();
private Map<DataTime, Collection<NPPSoundingRecord>> groupedRecords = new HashMap<DataTime, Collection<NPPSoundingRecord>>();
@ -113,11 +115,10 @@ public class NPPSoundingMapResource extends
public synchronized void addRecords(PluginDataObject... records) {
for (PluginDataObject record : records) {
if (record instanceof NPPSoundingRecord) {
if (allRecords.contains(record) == false) {
allRecords.add((NPPSoundingRecord) record);
}
}
}
Map<DataTime, Collection<NPPSoundingRecord>> groupedRecords = resourceData
.groupRecordTimes(allRecords);
List<DataTime> dataTimes = new ArrayList<DataTime>(
@ -184,7 +185,8 @@ public class NPPSoundingMapResource extends
*
*/
public boolean isEditable() {
return getCapability(EditableCapability.class).isEditable();
return getCapability(EditableCapability.class).isEditable()
&& getProperties().isVisible();
}
double getRadius() {

View file

@ -23,12 +23,11 @@ import java.text.DecimalFormat;
import java.text.ParsePosition;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.HashMap;
import java.util.Iterator;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
import javax.measure.Measure;
import javax.measure.converter.UnitConverter;
@ -36,7 +35,10 @@ import javax.measure.unit.Unit;
import javax.measure.unit.UnitFormat;
import org.geotools.coverage.grid.GeneralGridGeometry;
import org.opengis.geometry.Envelope;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.operation.MathTransform;
import com.raytheon.uf.common.colormap.prefs.ColorMapParameters;
import com.raytheon.uf.common.dataplugin.PluginDataObject;
@ -44,18 +46,18 @@ import com.raytheon.uf.common.dataplugin.npp.viirs.VIIRSDataRecord;
import com.raytheon.uf.common.datastorage.Request;
import com.raytheon.uf.common.datastorage.records.IDataRecord;
import com.raytheon.uf.common.datastorage.records.ShortDataRecord;
import com.raytheon.uf.common.geospatial.ISpatialObject;
import com.raytheon.uf.common.geospatial.ReferencedCoordinate;
import com.raytheon.uf.common.geospatial.util.EnvelopeIntersection;
import com.raytheon.uf.common.status.IUFStatusHandler;
import com.raytheon.uf.common.status.UFStatus;
import com.raytheon.uf.common.status.UFStatus.Priority;
import com.raytheon.uf.common.time.DataTime;
import com.raytheon.uf.common.time.TimeRange;
import com.raytheon.uf.viz.core.DrawableImage;
import com.raytheon.uf.viz.core.IGraphicsTarget;
import com.raytheon.uf.viz.core.IGraphicsTarget.RasterMode;
import com.raytheon.uf.viz.core.IMesh;
import com.raytheon.uf.viz.core.PixelCoverage;
import com.raytheon.uf.viz.core.PixelExtent;
import com.raytheon.uf.viz.core.datastructure.DataCubeContainer;
import com.raytheon.uf.viz.core.drawables.ColorMapLoader;
import com.raytheon.uf.viz.core.drawables.IImage;
@ -77,10 +79,16 @@ import com.raytheon.uf.viz.core.style.StyleRule;
import com.raytheon.uf.viz.core.tile.Tile;
import com.raytheon.uf.viz.core.tile.TileSetRenderable;
import com.raytheon.uf.viz.core.tile.TileSetRenderable.TileImageCreator;
import com.raytheon.uf.viz.npp.viirs.Activator;
import com.raytheon.uf.viz.npp.viirs.style.VIIRSDataRecordCriteria;
import com.raytheon.viz.core.style.image.DataScale;
import com.raytheon.viz.core.style.image.ImagePreferences;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Point;
import com.vividsolutions.jts.geom.prep.PreparedGeometry;
import com.vividsolutions.jts.geom.prep.PreparedGeometryFactory;
/**
* NPP VIIRS resource. Responsible for drawing a single color band.
@ -95,6 +103,7 @@ import com.vividsolutions.jts.geom.Coordinate;
* Nov 30, 2011 mschenke Initial creation
* Feb 21, 2012 #30 mschenke Fixed sampling issue
* Aug 2, 2013 #2190 mschenke Switched interrogate to use Measure objects
* Aug 27, 2013 #2190 mschenke Made interrogate more efficient
*
* </pre>
*
@ -142,17 +151,108 @@ public class VIIRSResource extends
/**
* Every {@link VIIRSDataRecord} will have a corresponding RecordData object
*/
private static class RecordData {
private class RecordData {
/** Indicates if data intersects current descriptor */
private boolean intersects = false;
private double INTERSECTION_FACTOR = 10.0;
/** Indicates if the {@link #tileSet} needs to be projected */
private boolean project = true;
/** Intersection geometry for the target */
private List<PreparedGeometry> targetIntersection;
/** Flag designated if a project call is required next paint */
private boolean project;
/** Renderable for the data record */
private TileSetRenderable tileSet;
private final double resolution;
public RecordData(VIIRSDataRecord dataRecord) {
this.tileSet = new TileSetRenderable(
getCapability(ImagingCapability.class), dataRecord
.getCoverage().getGridGeometry(),
new VIIRSTileImageCreator(dataRecord),
dataRecord.getLevels(), 256);
this.resolution = Math.min(dataRecord.getCoverage().getDx(),
dataRecord.getCoverage().getDy()) * INTERSECTION_FACTOR;
this.project = true;
}
public Collection<DrawableImage> getImagesToRender(
IGraphicsTarget target, PaintProperties paintProps)
throws VizException {
if (project) {
projectInternal();
project = false;
}
if (targetIntersection != null) {
return tileSet.getImagesToRender(target, paintProps);
} else {
return Collections.emptyList();
}
}
public void project() {
this.project = true;
}
private void projectInternal() {
GeneralGridGeometry targetGeometry = descriptor.getGridGeometry();
if (tileSet.getTargetGeometry() != targetGeometry) {
tileSet.project(targetGeometry);
try {
Envelope tileSetEnvelope = tileSet.getTileSetGeometry()
.getEnvelope();
targetIntersection = null;
Geometry intersection = EnvelopeIntersection
.createEnvelopeIntersection(
tileSetEnvelope,
targetGeometry.getEnvelope(),
resolution,
(int) (tileSetEnvelope.getSpan(0) / (resolution * INTERSECTION_FACTOR)),
(int) (tileSetEnvelope.getSpan(1) / (resolution * INTERSECTION_FACTOR)));
if (intersection != null) {
int numGeoms = intersection.getNumGeometries();
targetIntersection = new ArrayList<PreparedGeometry>(
numGeoms);
for (int n = 0; n < numGeoms; ++n) {
targetIntersection.add(PreparedGeometryFactory
.prepare(intersection.getGeometryN(n)
.buffer(resolution
* INTERSECTION_FACTOR)));
}
}
} catch (Exception e) {
Activator.statusHandler
.handle(Priority.PROBLEM,
"Error constructing envelope intersection for tileset",
e);
}
}
}
public double interrogate(Coordinate latLon) throws VizException {
return tileSet.interrogate(latLon);
}
public boolean contains(Geometry geom) {
if (targetIntersection != null) {
for (PreparedGeometry pg : targetIntersection) {
if (pg.contains(geom)) {
return true;
}
}
}
return false;
}
public void dispose() {
tileSet.dispose();
tileSet = null;
targetIntersection = null;
}
}
/**
@ -161,6 +261,8 @@ public class VIIRSResource extends
*/
private Map<VIIRSDataRecord, RecordData> dataRecordMap;
private Map<DataTime, Collection<VIIRSDataRecord>> groupedRecords;
private String name;
/**
@ -171,6 +273,7 @@ public class VIIRSResource extends
LoadProperties loadProperties, List<VIIRSDataRecord> initialRecords) {
super(resourceData, loadProperties);
dataRecordMap = new LinkedHashMap<VIIRSDataRecord, RecordData>();
groupedRecords = new HashMap<DataTime, Collection<VIIRSDataRecord>>();
resourceData.addChangeListener(this);
dataTimes = new ArrayList<DataTime>();
for (VIIRSDataRecord record : initialRecords) {
@ -208,25 +311,8 @@ public class VIIRSResource extends
RecordData data = dataRecordMap.get(dataRecord);
if (data == null) {
// New record
data = new RecordData();
data.tileSet = new TileSetRenderable(
getCapability(ImagingCapability.class), dataRecord
.getCoverage().getGridGeometry(),
new VIIRSTileImageCreator(dataRecord),
dataRecord.getLevels(), 256);
dataRecordMap.put(dataRecord, data);
dataRecordMap.put(dataRecord, new RecordData(dataRecord));
}
// Make sure record intersects map
data.intersects = intersects(dataRecord);
}
private boolean intersects(VIIRSDataRecord dataRecord) {
// Make sure spatial record intersects map
PixelCoverage dataCoverage = descriptor.worldToPixel(dataRecord
.getCoverage().getEnvelope().getEnvelopeInternal());
return dataCoverage.intersects(new PixelExtent(descriptor
.getGridGeometry().getGridRange()));
}
/**
@ -374,6 +460,32 @@ public class VIIRSResource extends
colorMapParameters);
}
public void addRecords(PluginDataObject... records) {
synchronized (dataRecordMap) {
for (PluginDataObject record : records) {
if (record instanceof VIIRSDataRecord) {
try {
addRecord((VIIRSDataRecord) record);
} catch (VizException e) {
statusHandler.handle(
Priority.PROBLEM,
"Error adding record from update: "
+ e.getLocalizedMessage(), e);
}
}
}
Map<DataTime, Collection<VIIRSDataRecord>> groupedRecords = resourceData
.groupRecordTimes(dataRecordMap.keySet());
List<DataTime> dataTimes = new ArrayList<DataTime>(
groupedRecords.keySet());
Collections.sort(dataTimes);
this.dataTimes = dataTimes;
this.groupedRecords = groupedRecords;
}
}
/*
* (non-Javadoc)
*
@ -383,22 +495,20 @@ public class VIIRSResource extends
*/
@Override
public void remove(DataTime dataTime) {
super.remove(dataTime);
synchronized (dataRecordMap) {
Iterator<Entry<VIIRSDataRecord, RecordData>> iter = dataRecordMap
.entrySet().iterator();
while (iter.hasNext()) {
Entry<VIIRSDataRecord, RecordData> entry = iter.next();
VIIRSDataRecord record = entry.getKey();
if (VIIRSResourceData.withinRange(dataTime.getValidPeriod(),
record.getDataTime())) {
iter.remove();
RecordData data = entry.getValue();
data.tileSet.dispose();
Collection<VIIRSDataRecord> records = groupedRecords
.remove(dataTime);
if (records != null) {
for (VIIRSDataRecord record : records) {
RecordData data = dataRecordMap.remove(record);
if (data != null) {
data.dispose();
}
}
}
}
super.remove(dataTime);
}
/*
* (non-Javadoc)
@ -409,9 +519,10 @@ public class VIIRSResource extends
protected void disposeInternal() {
synchronized (dataRecordMap) {
for (RecordData data : dataRecordMap.values()) {
data.tileSet.dispose();
data.dispose();
}
dataRecordMap.clear();
groupedRecords.clear();
dataTimes.clear();
}
}
@ -447,8 +558,9 @@ public class VIIRSResource extends
protected void initInternal(IGraphicsTarget target) throws VizException {
getCapability(ImagingCapability.class).setProvider(this);
synchronized (dataRecordMap) {
resourceChanged(ChangeType.DATA_UPDATE, dataRecordMap.keySet()
.toArray(new PluginDataObject[dataRecordMap.size()]));
PluginDataObject[] pdos = dataRecordMap.keySet().toArray(
new PluginDataObject[0]);
resourceChanged(ChangeType.DATA_UPDATE, pdos);
}
}
@ -468,25 +580,8 @@ public class VIIRSResource extends
for (int i = 0; i < pdos.length; ++i) {
records[i] = (VIIRSDataRecord) pdos[i];
}
synchronized (dataRecordMap) {
long t0 = System.currentTimeMillis();
for (VIIRSDataRecord record : records) {
try {
addRecord(record);
} catch (VizException e) {
statusHandler.handle(
Priority.PROBLEM,
"Error adding record from update: "
+ e.getLocalizedMessage(), e);
}
}
System.out.println("Time to add " + records.length
+ " records: " + (System.currentTimeMillis() - t0)
+ "ms");
dataTimes = new ArrayList<DataTime>(resourceData
.groupRecordTimes(dataRecordMap.keySet()).keySet());
}
addRecords(records);
}
issueRefresh();
}
@ -537,38 +632,54 @@ public class VIIRSResource extends
Unit<?> dataUnit = params.getDisplayUnit();
double noDataValue = params.getNoDataValue();
Coordinate latLon = null;
Coordinate crs = null;
try {
latLon = coord.asLatLon();
Coordinate grid = coord.asGridCell(descriptor.getGridGeometry(),
PixelInCell.CELL_CENTER);
MathTransform mt = descriptor.getGridGeometry().getGridToCRS(
PixelInCell.CELL_CENTER);
double[] out = new double[2];
mt.transform(new double[] { grid.x, grid.y }, 0, out, 0, 1);
crs = new Coordinate(out[0], out[1]);
} catch (Exception e) {
throw new VizException(
"Could not get lat/lon from ReferencedCoordinate", e);
}
DataTime currTime = descriptor.getTimeForResource(this);
if (currTime != null) {
Point crsPoint = new GeometryFactory().createPoint(crs);
double bestValue = Double.NaN;
// Since there is overlap between granules, the best value is the
// one that is closest to 0
VIIRSDataRecord bestRecord = null;
synchronized (dataRecordMap) {
for (VIIRSDataRecord record : dataRecordMap.keySet()) {
if (VIIRSResourceData.withinRange(
currTime.getValidPeriod(), record.getDataTime())) {
Collection<VIIRSDataRecord> records = groupedRecords.get(descriptor
.getTimeForResource(this));
if (records != null) {
// Since there is overlap between granules, the best value is
// the one that is closest to 0
for (VIIRSDataRecord record : records) {
RecordData data = dataRecordMap.get(record);
if (data.intersects && data.project == false) {
double value = data.tileSet.interrogate(latLon);
if (data.contains(crsPoint)) {
double value = data.interrogate(latLon);
if (Double.isNaN(value) == false
&& value != noDataValue) {
bestValue = value;
bestRecord = record;
break;
}
}
}
}
}
if (Double.isNaN(bestValue) == false) {
dataValue = params.getDataToDisplayConverter().convert(
bestValue);
}
if (bestRecord != null) {
interMap.put(ISpatialObject.class.toString(),
bestRecord.getSpatialObject());
}
}
interMap.put(SATELLITE_DATA_INTERROGATE_ID,
@ -587,10 +698,8 @@ public class VIIRSResource extends
public void project(CoordinateReferenceSystem crs) throws VizException {
synchronized (dataRecordMap) {
try {
for (VIIRSDataRecord record : dataRecordMap.keySet()) {
RecordData data = dataRecordMap.get(record);
data.intersects = intersects(record);
data.project = true;
for (RecordData data : dataRecordMap.values()) {
data.project();
}
} catch (Exception e) {
statusHandler.handle(Priority.PROBLEM,
@ -611,22 +720,15 @@ public class VIIRSResource extends
public Collection<DrawableImage> getImages(IGraphicsTarget target,
PaintProperties paintProps) throws VizException {
List<DrawableImage> images = new ArrayList<DrawableImage>();
DataTime currTime = paintProps.getDataTime();
if (currTime != null) {
TimeRange tr = currTime.getValidPeriod();
synchronized (dataRecordMap) {
for (VIIRSDataRecord record : dataRecordMap.keySet()) {
if (VIIRSResourceData.withinRange(tr, record.getDataTime())) {
Collection<VIIRSDataRecord> records = groupedRecords.get(paintProps
.getDataTime());
if (records != null) {
for (VIIRSDataRecord record : records) {
RecordData data = dataRecordMap.get(record);
if (data.intersects) {
if (data.project) {
data.tileSet.project(descriptor
.getGridGeometry());
data.project = false;
}
images.addAll(data.tileSet.getImagesToRender(
target, paintProps));
}
if (data != null) {
images.addAll(data
.getImagesToRender(target, paintProps));
}
}
}

View file

@ -19,16 +19,15 @@
**/
package com.raytheon.uf.viz.npp;
import java.util.ArrayDeque;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
import java.util.HashMap;
import java.util.Iterator;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Queue;
import java.util.Set;
import javax.xml.bind.annotation.XmlAccessType;
import javax.xml.bind.annotation.XmlAccessorType;
@ -39,6 +38,7 @@ import com.raytheon.uf.common.dataquery.requests.RequestConstraint;
import com.raytheon.uf.common.dataquery.requests.RequestConstraint.ConstraintType;
import com.raytheon.uf.common.time.DataTime;
import com.raytheon.uf.common.time.TimeRange;
import com.raytheon.uf.common.time.util.TimeUtil;
import com.raytheon.uf.viz.core.catalog.LayerProperty;
import com.raytheon.uf.viz.core.datastructure.DataCubeContainer;
import com.raytheon.uf.viz.core.exception.VizException;
@ -56,6 +56,7 @@ import com.raytheon.uf.viz.core.rsc.AbstractRequestableResourceData;
* ------------ ---------- ----------- --------------------------
* Jan 8, 2013 mschenke Initial creation
* Aug 2, 2013 2190 mschenke Moved time grouping functions to utility class
* Aug 27, 2013 2190 mschenke Fixed groupRecordTimes function
*
* </pre>
*
@ -66,34 +67,6 @@ import com.raytheon.uf.viz.core.rsc.AbstractRequestableResourceData;
public abstract class AbstractNppResourceData extends
AbstractRequestableResourceData {
private static class DataTimeIterator<T extends PluginDataObject>
implements Iterator<DataTime> {
private Iterator<T> pdoIter;
private T lastAccesed;
private DataTimeIterator(Iterator<T> pdoIter) {
this.pdoIter = pdoIter;
}
@Override
public boolean hasNext() {
return pdoIter.hasNext();
}
@Override
public DataTime next() {
lastAccesed = pdoIter.next();
return lastAccesed.getDataTime();
}
@Override
public void remove() {
pdoIter.remove();
}
}
@XmlElement
private int groupTimeRangeMinutes = 15;
@ -129,14 +102,17 @@ public abstract class AbstractNppResourceData extends
getMetadataMap());
RequestConstraint timeConst = new RequestConstraint();
timeConst.setConstraintType(ConstraintType.BETWEEN);
timeConst.setBetweenValueList(new String[] {
new DataTime(first.getValidPeriod().getStart()).toString(),
new DataTime(last.getValidPeriod().getEnd()).toString() });
timeConst
.setBetweenValueList(new String[] {
TimeUtil.formatToSqlTimestamp(first.getValidPeriod()
.getStart()),
TimeUtil.formatToSqlTimestamp(last.getValidPeriod()
.getEnd()) });
requestMap.put("dataTime.refTime", timeConst);
LayerProperty property = new LayerProperty();
property.setEntryQueryParameters(requestMap, false);
property.setNumberOfImages(9999);
property.setNumberOfImages(Integer.MAX_VALUE);
List<Object> pdos = DataCubeContainer.getData(property, 60000);
List<PluginDataObject> finalList = new ArrayList<PluginDataObject>(
@ -152,8 +128,8 @@ public abstract class AbstractNppResourceData extends
}
}
}
Collections.sort(finalList, layerComparator);
}
Collections.sort(finalList, layerComparator);
return finalList.toArray(new PluginDataObject[finalList.size()]);
}
@ -185,29 +161,26 @@ public abstract class AbstractNppResourceData extends
*/
public <T extends PluginDataObject> Map<DataTime, Collection<T>> groupRecordTimes(
Collection<T> records) {
long groupTimeInMillis = groupTimeRangeMinutes * 60 * 1000;
Map<DataTime, Collection<T>> grouped = new HashMap<DataTime, Collection<T>>();
Queue<T> objects = new ArrayDeque<T>(records);
while (objects.size() > 0) {
T record = objects.remove();
DataTime current = record.getDataTime();
TimeRange prev, curr;
prev = curr = current.getValidPeriod();
List<T> group = new ArrayList<T>();
while (curr != null) {
prev = curr;
DataTimeIterator<T> iter = new DataTimeIterator<T>(
objects.iterator());
curr = NPPTimeUtility.match(iter, prev, groupTimeInMillis);
if (curr != null) {
group.add(iter.lastAccesed);
}
Set<DataTime> times = new HashSet<DataTime>();
for (T record : records) {
times.add(record.getDataTime());
}
grouped.put(new DataTime(prev.getStart().getTime(), prev), group);
Collection<DataTime> groupedTimes = groupTimes(times);
Map<DataTime, List<T>> grouped = new HashMap<DataTime, List<T>>();
for (T record : records) {
for (DataTime time : groupedTimes) {
if (withinRange(time.getValidPeriod(), record.getDataTime())) {
List<T> group = grouped.get(time);
if (group == null) {
group = new ArrayList<T>();
grouped.put(time, group);
}
return grouped;
group.add(record);
}
}
}
return new HashMap<DataTime, Collection<T>>(grouped);
}
/**

View file

@ -102,6 +102,8 @@ public class PopupSkewTContainer {
this.soundingProvider = VerticalSoundingProviderFactory
.getVerticalSoundingProvider(soundingSource.getType(),
soundingSource.getConstraints());
} else {
this.soundingProvider = null;
}
}

View file

@ -273,7 +273,7 @@ public class PlotModelDataRequestJob extends Job {
if (pdc == null) {
// Datacube didn't have proper plugin; going
// directly to the data store
pdc = PointDataRequest.requestPointDataAllLevels(null,
pdc = PointDataRequest.requestPointDataAllLevels(
this.plugin,
params.toArray(new String[params.size()]), null,
map);

View file

@ -31,6 +31,7 @@ import com.raytheon.uf.common.pointdata.PointDataContainer;
import com.raytheon.uf.common.pointdata.PointDataServerRequest;
import com.raytheon.uf.common.time.DataTime;
import com.raytheon.uf.common.time.TimeRange;
import com.raytheon.uf.common.time.util.TimeUtil;
import com.raytheon.uf.viz.core.exception.VizException;
import com.raytheon.uf.viz.core.requests.ThriftClient;
@ -63,6 +64,8 @@ public class PointDataRequest {
private static final String DATATIME_KEY = "dataTime";
private static final String REFTIME_KEY = DATATIME_KEY + ".refTime";
private static final String REQUESTED_PARAMETERS_KEY = "requestedParameters";
private static final String RESTRICT_LEVEL = "restrictLevel";
@ -120,6 +123,30 @@ public class PointDataRequest {
stationIds, constraints, true, null, null);
}
/**
* Request all levels for point data for a set of stations or points, over a
* given set of parameters
*
* The request can be additionally constrained by optional
* RequestConstraints.
*
* @param pluginName
* the plugin to use (required)
* @param parameters
* the parameters to request (required)
* @param stationIds
* the station IDs to constrain to (optional)
* @param constraints
* additional constraints (optional)
* @return
*/
public static PointDataContainer requestPointDataAllLevels(
String pluginName, String[] parameters, String[] stationIds,
Map<String, RequestConstraint> constraints) throws VizException {
return requestPointDataInternal(null, null, pluginName, parameters,
stationIds, constraints, true, null, null);
}
/**
* Request all levels for point data for a set of stations or points, over a
* given set of parameters, at a specific point in time.
@ -150,6 +177,32 @@ public class PointDataRequest {
stationIds, constraints, true, null, null);
}
/**
* Request all levels for point data for a set of stations or points, over a
* given set of parameters, over a time range.
*
* The request can be additionally constrained by optional
* RequestConstraints.
*
* @param tr
* the time range to request the data for (required)
* @param pluginName
* the plugin to use (required)
* @param parameters
* the parameters to request (required)
* @param stationIds
* the station IDs to constrain to (optional)
* @param constraints
* additional constraints (optional)
* @return
*/
public static PointDataContainer requestPointDataAllLevels(TimeRange tr,
String pluginName, String[] parameters, String[] stationIds,
Map<String, RequestConstraint> constraints) throws VizException {
return requestPointDataInternal(null, tr, pluginName, parameters,
stationIds, constraints, true, null, null);
}
public static String[] getParameterNames(String pluginName,
Map<String, RequestConstraint> constraints) throws VizException {
try {
@ -240,14 +293,13 @@ public class PointDataRequest {
dtConstraint.setConstraintType(ConstraintType.IN);
rcMap.put(DATATIME_KEY, dtConstraint);
} else if (tr != null) {
DataTime start = new DataTime(tr.getStart());
DataTime end = new DataTime(tr.getEnd());
RequestConstraint dtConstraint = new RequestConstraint();
String[] constraintList = { start.toString(), end.toString() };
String[] constraintList = {
TimeUtil.formatToSqlTimestamp(tr.getStart()),
TimeUtil.formatToSqlTimestamp(tr.getEnd()) };
dtConstraint.setBetweenValueList(constraintList);
dtConstraint.setConstraintType(ConstraintType.BETWEEN);
rcMap.put(DATATIME_KEY, dtConstraint);
rcMap.put(REFTIME_KEY, dtConstraint);
}
if (stationIds != null) {

View file

@ -235,7 +235,7 @@ public class PointDataCubeAdapter implements IDataCubeAdapter {
public PointDataContainer getBaseRecords(Collection<String> baseParams,
Map<String, RequestConstraint> queryParams) throws VizException {
String plugin = queryParams.get(PLUGIN_NAME).getConstraintValue();
return PointDataRequest.requestPointDataAllLevels(null, plugin,
return PointDataRequest.requestPointDataAllLevels(plugin,
baseParams.toArray(new String[] {}), null, queryParams);
}

View file

@ -20,7 +20,6 @@
package com.raytheon.uf.common.dataplugin.npp.viirs.projection;
import static java.lang.Math.abs;
import static java.lang.Math.acos;
import static java.lang.Math.atan;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
@ -33,7 +32,6 @@ import java.util.Arrays;
import org.geotools.parameter.DefaultParameterDescriptor;
import org.geotools.parameter.DefaultParameterDescriptorGroup;
import org.geotools.referencing.GeodeticCalculator;
import org.geotools.referencing.operation.MathTransformProvider;
import org.geotools.referencing.operation.projection.MapProjection;
import org.geotools.referencing.operation.projection.ProjectionException;
@ -45,13 +43,12 @@ import org.opengis.parameter.ParameterNotFoundException;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.LineSegment;
/**
* TODO Add Description
* VIIRS map projection
*
* <pre>
*
@ -64,6 +61,7 @@ import com.vividsolutions.jts.geom.LineSegment;
* calculations. Switched projection extrapolation
* algorithm to use great circle intersections for
* more accurate results
* Aug 27, 2013 #2190 mschenke Sped up transform functions
*
* </pre>
*
@ -107,15 +105,17 @@ public class VIIRSMapProjection extends MapProjection {
public static final double PIO2 = Math.PI / 2;
private float[] centerLats;
private final float[] centerLats;
private float[] centerLons;
private final float[] centerLons;
private float[] directions;
private final float[] directions;
private double resolution;
private final double resolution;
private int actualHeight;
private final int validHeight;
private final int actualHeight;
/**
* @param values
@ -129,6 +129,12 @@ public class VIIRSMapProjection extends MapProjection {
this.directions = Provider.getValue(Provider.CENTER_DIRECTIONS, values);
this.resolution = Provider.getValue(Provider.RESOLUTION, values);
this.actualHeight = Provider.getValue(Provider.CENTER_LENGTH, values);
if (centerLats.length != centerLons.length
|| centerLats.length != directions.length) {
throw new IllegalArgumentException(
"Center lat/lon and direction arrays must be same length");
}
this.validHeight = centerLats.length;
}
/*
@ -153,115 +159,145 @@ public class VIIRSMapProjection extends MapProjection {
Point2D ptDst) throws ProjectionException {
double xi = x / resolution;
double yi = y / resolution;
double lon = Math.toRadians(getInterpolatedValue(xi - 0.5, actualHeight
- yi - 0.5, 0));
double lat = Math.toRadians(getInterpolatedValue(xi - 0.5, actualHeight
- yi - 0.5, 1));
Coordinate latLon = getInterpolatedCoordinate(xi - 0.5, actualHeight
- yi - 0.5);
Point2D point = ptDst != null ? ptDst : new Point2D.Double();
point.setLocation(lon, lat);
point.setLocation(Math.toRadians(latLon.x), Math.toRadians(latLon.y));
return point;
}
protected float getInterpolatedValue(double xd, double yd, int dim) {
float baseVal = Float.NaN;
float value = 0.0f;
float missing = 1.0f;
/**
* @param d
* @param e
* @return
*/
private Coordinate getInterpolatedCoordinate(double xd, double yd) {
double baseLon = Double.NaN;
float latValue = 0.0f;
float lonValue = 0.0f;
float latMissing = 1.0f;
float lonMissing = 1.0f;
float x = (float) xd;
float y = (float) yd;
int xi = (int) x;
int yi = (int) y;
int xi = (int) xd;
int yi = (int) yd;
// Upper left corner
float xWeight = 1 - x + xi;
float yWeight = 1 - y + yi;
float weight = xWeight * yWeight;
float val = getRawDataValue(xi, yi, dim);
if (Float.isNaN(val)) {
missing -= weight;
Coordinate c = getRawDataCoordinate(xi, yi);
if (Double.isNaN(c.x)) {
lonMissing -= weight;
} else {
value += weight * val;
baseVal = val;
lonValue += weight * c.x;
baseLon = c.x;
}
if (Double.isNaN(c.y)) {
latMissing -= weight;
} else {
latValue += weight * c.y;
}
// upper right corner
xi = xi + 1;
xWeight = 1 - xi + x;
yWeight = 1 - y + yi;
weight = xWeight * yWeight;
val = getRawDataValue(xi, yi, dim);
if (Float.isNaN(val)) {
missing -= weight;
c = getRawDataCoordinate(xi, yi);
if (Double.isNaN(c.x)) {
lonMissing -= weight;
} else {
if (Float.isNaN(baseVal)) {
baseVal = val;
if (Double.isNaN(baseLon)) {
baseLon = c.x;
} else {
val = correct(baseVal, val, dim);
c.x = correct(baseLon, c.x);
}
value += weight * val;
lonValue += weight * c.x;
}
if (Double.isNaN(c.y)) {
latMissing -= weight;
} else {
latValue += weight * c.y;
}
// lower right corner
yi = yi + 1;
xWeight = 1 - xi + x;
yWeight = 1 - yi + y;
weight = xWeight * yWeight;
val = getRawDataValue(xi, yi, dim);
if (Float.isNaN(val)) {
missing -= weight;
c = getRawDataCoordinate(xi, yi);
if (Double.isNaN(c.x)) {
lonMissing -= weight;
} else {
if (Float.isNaN(baseVal)) {
baseVal = val;
if (Double.isNaN(baseLon)) {
baseLon = c.x;
} else {
val = correct(baseVal, val, dim);
c.x = correct(baseLon, c.x);
}
value += weight * val;
lonValue += weight * c.x;
}
if (Double.isNaN(c.y)) {
latMissing -= weight;
} else {
latValue += weight * c.y;
}
// lower left corner
xi = xi - 1;
xWeight = 1 - x + xi;
yWeight = 1 - yi + y;
weight = xWeight * yWeight;
val = getRawDataValue(xi, yi, dim);
if (Float.isNaN(val)) {
missing -= weight;
c = getRawDataCoordinate(xi, yi);
if (Double.isNaN(c.x)) {
lonMissing -= weight;
} else {
if (Float.isNaN(baseVal) == false) {
val = correct(baseVal, val, dim);
if (Double.isNaN(baseLon) == false) {
c.x = correct(baseLon, c.x);
}
value += weight * val;
lonValue += weight * c.x;
}
if (Double.isNaN(c.y)) {
latMissing -= weight;
} else {
latValue += weight * c.y;
}
value /= missing;
if (dim == 0) {
lonValue /= lonMissing;
latValue /= latMissing;
// If we did any normalization, reverse it here
if (value > 180.0f) {
value -= 360.0f;
} else if (value <= -180.0f) {
value += 360.0f;
}
}
return value / missing;
if (lonValue > 180.0f) {
lonValue -= 360.0f;
} else if (lonValue <= -180.0f) {
lonValue += 360.0f;
}
private float correct(double baseVal, double val, int dim) {
if (dim == 0 && Math.abs(baseVal - val) > 180.0f) {
if (baseVal > val) {
val += 360.0f;
return new Coordinate(lonValue, latValue);
}
private double correct(double baseLon, double lon) {
if (Math.abs(baseLon - lon) > 180.0f) {
if (baseLon > lon) {
lon += 360.0f;
} else {
val -= 360.0f;
lon -= 360.0f;
}
}
return (float) val;
return lon;
}
/**
* @param xi
* @param yi
* @param dim
* @return
*/
private float getRawDataValue(int xi, int yi, int dim) {
private Coordinate getRawDataCoordinate(int xi, int yi) {
Coordinate c = null;
if (yi >= 0 && yi < centerLats.length) {
if (yi >= 0 && yi < validHeight) {
c = index(xi, yi);
} else {
int closestY = yi;
@ -270,7 +306,7 @@ public class VIIRSMapProjection extends MapProjection {
closestY = 0;
closestY2 = closestY + 1;
} else {
closestY = centerLats.length - 1;
closestY = validHeight - 1;
closestY2 = closestY - 1;
}
@ -280,9 +316,8 @@ public class VIIRSMapProjection extends MapProjection {
Coordinate b = index(xi, closestY2);
LineSegment ls = new LineSegment(a, b);
c = ls.pointAlong(-Math.abs(yDiff));
}
return (float) (dim == 0 ? c.x : c.y);
return c;
}
private Coordinate index(int xi, int yi) {
@ -298,7 +333,154 @@ public class VIIRSMapProjection extends MapProjection {
double earthRadius = getRadius(latitude);
double[] latLon = target_pt(latitude, Math.toRadians(nadirCoord.x),
(resolution * abs(xi)) / earthRadius, az);
return new Coordinate(latLon[0], latLon[1]);
return new Coordinate(Math.toDegrees(latLon[0]),
Math.toDegrees(latLon[1]));
}
private static class DistanceResult {
private final int index;
private final double xDist;
private final double yDist;
public DistanceResult(int index, double xDist, double yDist) {
this.index = index;
this.xDist = xDist;
this.yDist = yDist;
}
}
private DistanceResult distance(int index, double aLonR, double aLatR) {
double clon = centerLons[index];
double clat = centerLats[index];
// Compute azm_sidb for center point to lat/lon point
double[] azm_sidb = azm_sidb(Math.toRadians(clat),
Math.toRadians(clon), aLatR, aLonR);
double actual = azm_sidb[0];
double expected = directions[index];
// Correct for world wrapping
if (Math.abs(actual - expected) > Math.PI) {
if (actual > expected) {
expected += TWOPI;
} else {
expected -= TWOPI;
}
}
double xFactor = 1;
double yFactor = 1;
if (actual < expected) {
expected -= PIO2;
yFactor = -1;
} else {
expected += PIO2;
xFactor = -1;
}
// Get angle difference between actual and expected as well as
// multiplication factor for whether left or right of center points
double diff = actual - expected;
// Because diff is always [0,90], we must account for what quadrant it
// is in and adjust sign appropriately based on factor
double yDist = yFactor * azm_sidb[1] * Math.sin(diff);
double xDist = xFactor * azm_sidb[1] * Math.cos(diff);
return new DistanceResult(index, xDist, yDist);
}
/**
* Intelligently searches for the best index for the given lat/lon in the
* range of startIdx->endIdx
*
* @param startIdx
* @param endIdx
* @param aLonR
* @param aLatR
* @return
*/
private DistanceResult findBest(int startIdx, int endIdx, double aLonR,
double aLatR) {
// Get search results for start/end idx
DistanceResult below = distance(startIdx, aLonR, aLatR);
DistanceResult above = distance(endIdx, aLonR, aLatR);
DistanceResult best = null;
boolean done = false;
do {
double aboveDist = above.yDist;
double belowDist = below.yDist;
int aboveIndex = above.index;
int belowIndex = below.index;
// get range of y distances
double range = aboveDist - belowDist;
// Calculate best index to search when looking for yDist of 0
double testIdx = belowIndex + (-belowDist / range)
* (aboveIndex - belowIndex);
// Ensure index within search bounds
int index = (int) testIdx;
if (index < belowIndex) {
index = belowIndex;
} else if (index > aboveIndex) {
index = aboveIndex;
}
// Get other index to check which will be 1 greater than index due
// to casting of testIdx
int otherIndex = index + 1;
if (otherIndex >= validHeight) {
otherIndex = validHeight - 1;
}
if (index == belowIndex || index == aboveIndex
|| otherIndex == aboveIndex) {
// Exit case, one of the indices to search is same as our bounds
if (index == belowIndex) {
best = below;
} else if (index == aboveIndex) {
best = above;
} else {
// otherIndex is equal to above index, swap index and
// otherIndex and search index before other
index = otherIndex;
otherIndex = index - 1;
best = above;
}
// Search other index, set best result
DistanceResult other = distance(otherIndex, aLonR, aLatR);
double otherDist = Math.abs(other.yDist);
if (otherDist < Math.abs(best.yDist)) {
best = other;
}
// Mark done
done = true;
} else {
// Search the test index and determine new above/below indices
DistanceResult test = distance(index, aLonR, aLatR);
double testDist = test.yDist;
if (testDist < 0) {
if (aboveDist > 0) {
above = test;
} else {
below = test;
}
} else if (testDist > 0) {
if (aboveDist < 0) {
below = test;
} else {
above = test;
}
} else {
// Exactly 0!!! we are done!!
best = test;
done = true;
}
}
} while (!done);
return best;
}
/*
@ -311,371 +493,24 @@ public class VIIRSMapProjection extends MapProjection {
protected Point2D transformNormalized(double aLonR, double aLatR,
Point2D ptDst) throws ProjectionException {
double bestY = Double.NaN, bestX = 1;
DistanceResult best = findBest(0, validHeight - 1, aLonR, aLatR);
int idxToUse = best.index;
// Get radius to convert xDist/yDist into meters
double radius = getRadius(Math.toRadians(centerLats[idxToUse]));
Coordinate c1 = new Coordinate(Math.toRadians(centerLons[0]),
Math.toRadians(centerLats[0]));
Coordinate c2 = new Coordinate(Math.toRadians(centerLons[1]),
Math.toRadians(centerLats[1]));
// Compute bestX value based on hypotenuse and angle diff
bestX = resolution * 0.5 + best.xDist * radius;
// Extrapolate from the bottom of the image
double[] bottom = extrapolate(c1, c2, aLonR, aLatR);
// Compute potential bestY grid value
bestY = actualHeight - idxToUse - 0.5;
// Compute bestY value based on hypotenuse and angle diff in grid space
bestY = bestY * resolution - best.yDist * radius;
c1 = new Coordinate(Math.toRadians(centerLons[centerLons.length - 2]),
Math.toRadians(centerLats[centerLats.length - 2]));
c2 = new Coordinate(Math.toRadians(centerLons[centerLons.length - 1]),
Math.toRadians(centerLats[centerLats.length - 1]));
// Extrapolate from the top of the image
double[] top = extrapolate(c1, c2, aLonR, aLatR);
double dist1 = bottom[0], potentialY1;
double angleDiff1 = bottom[1];
double potentialX1 = bottom[2];
// Check to see if intersection point used is going same direction
// from c1 to c2 or if opposite
if (angleDiff1 > 90.0) {
potentialY1 = actualHeight + dist1;
} else {
potentialY1 = actualHeight - dist1;
}
double dist2 = top[0], potentialY2;
double angleDiff2 = top[1];
double potentialX2 = top[2];
// Check to see if intersection point used is going same direction
// from c1 to c2 or if opposite
if (angleDiff2 > 90.0) {
potentialY2 = actualHeight + dist2 - centerLons.length - 2;
} else {
potentialY2 = actualHeight - dist2 - centerLons.length - 2;
}
int validHeight = centerLats.length;
if ((potentialY1 >= 0 && potentialY1 < validHeight)
|| (potentialY2 >= 0 && potentialY2 < validHeight)) {
potentialY1 = validHeight - potentialY1;
potentialY2 = validHeight - potentialY2;
int startI = (int) Math.floor(Math.min(potentialY1, potentialY2));
if (startI < 0) {
startI = 0;
}
int endI = (int) Math.ceil(Math.max(potentialY1, potentialY2));
if (endI >= validHeight) {
endI = validHeight - 1;
}
// Lat/Lon to CRS, within our known area
double alat = Math.toDegrees(aLatR);
double alon = Math.toDegrees(aLonR);
double bestDiff = Double.NaN;
double bestDist = Double.NaN;
double bestFactor = 1;
for (int i = startI; i <= endI; ++i) {
double clon = centerLons[i];
double clat = centerLats[i];
double[] azm_sidb = azm_sidb(Math.toRadians(clat),
Math.toRadians(clon), aLatR, aLonR);
double actual = azm_sidb[0];
double expected = directions[i];
// Correct for world wrapping
if (Math.abs(actual - expected) > Math.PI) {
if (actual > expected) {
expected += TWOPI;
} else {
expected -= TWOPI;
}
}
double factor = 1;
if (actual < expected) {
expected -= PIO2;
} else {
expected += PIO2;
factor = -1;
}
double diff = Math.abs(actual - expected);
double dist = azm_sidb[1] * getRadius(Math.toRadians(clat));
if (Double.isNaN(bestDiff) || diff <= bestDiff) {
bestFactor = factor;
bestDist = dist;
bestDiff = diff;
bestY = actualHeight - i - 0.5;
}
}
// We found a point which meets our threshold criteria, this means
// that the point is somewhere between our valid Y range
bestDist = Math.cos(bestDiff) * bestDist;
bestX = 0.5 + (bestFactor * (bestDist / resolution));
int yIdx = (int) bestY;
int otherYIdx;
if (yIdx < centerLats.length - 1) {
otherYIdx = yIdx + 1;
} else {
otherYIdx = yIdx - 1;
}
try {
Point2D latLon = new Point2D.Double(), otherLatLon = new Point2D.Double();
inverse().transform(
new Point2D.Double(bestX * resolution, (yIdx + 0.5)
* resolution), latLon);
inverse().transform(
new Point2D.Double(bestX * resolution,
(otherYIdx + 0.5) * resolution), otherLatLon);
Coordinate a = new Coordinate(latLon.getX(), latLon.getY());
Coordinate b = new Coordinate(otherLatLon.getX(),
otherLatLon.getY());
// Correct for world wrapping a, b and alon
if (Math.abs(a.x - b.x) > 180.0) {
if (a.x > b.x) {
b.x += 360.0;
} else {
b.x -= 360.0;
}
}
if (Math.abs(a.x - alon) > 180.0) {
if (a.x > alon) {
alon += 360.0;
} else {
alon -= 360.0;
}
}
LineSegment ls = new LineSegment(a, b);
if (yIdx < centerLats.length - 1) {
bestY += ls.projectionFactor(new Coordinate(alon, alat));
} else {
bestY -= ls.projectionFactor(new Coordinate(alon, alat));
}
} catch (TransformException e) {
return null;
}
} else {
// Outside range, use closest distance for values
if (dist1 < dist2) {
bestX = potentialX1;
bestY = potentialY1;
} else {
bestX = potentialX2;
bestY = potentialY2;
}
}
Point2D point = ptDst != null ? ptDst : new Point2D.Double();
point.setLocation(bestX * resolution, bestY * resolution);
point.setLocation(bestX, bestY);
return point;
}
/**
* Given two coordinates (radians), intersect the great circle for those two
* coordinates with the great circle perpendicular that goes through point
* (aLonR, aLatR)
*
* @param c1
* @param c2
* @param aLonR
* @param aLatR
* @return the minimum dist value to c1 [0], the angle difference of the
* intersecting point and c1-c2 angle [1] and the x value [2]
*/
private double[] extrapolate(Coordinate c1, Coordinate c2, double aLonR,
double aLatR) {
// Create GeodeticCalculator for calculations...
GeodeticCalculator gc = new GeodeticCalculator();
// Set start/end points as c1/c2
gc.setStartingGeographicPoint(Math.toDegrees(c1.x),
Math.toDegrees(c1.y));
gc.setDestinationGeographicPoint(Math.toDegrees(c2.x),
Math.toDegrees(c2.y));
// Get the angle of
double c12Dir = gc.getAzimuth();
double c12Dist = gc.getOrthodromicDistance();
// Set starting point (c1) and calculate another point on the great
// circle which we will call (c3)
gc.setStartingGeographicPoint(Math.toDegrees(c1.x),
Math.toDegrees(c1.y));
gc.setDirection(c12Dir, 2 * c12Dist);
Point2D dest = gc.getDestinationGeographicPoint();
Coordinate c3 = new Coordinate(Math.toRadians(dest.getX()),
Math.toRadians(dest.getY()));
// Get angle from c2 to c3
gc.setStartingGeographicPoint(Math.toDegrees(c2.x),
Math.toDegrees(c2.y));
gc.setDestinationGeographicPoint(Math.toDegrees(c3.x),
Math.toDegrees(c3.y));
double c23Dir = gc.getAzimuth();
double c23Dist = gc.getOrthodromicDistance();
// Get point perpendicular to c1 (c1_90)
gc.setStartingGeographicPoint(Math.toDegrees(c1.x),
Math.toDegrees(c1.y));
gc.setDirection(adjustAngle(c12Dir + 90.0), c12Dist);
dest = gc.getDestinationGeographicPoint();
Coordinate c1_90 = new Coordinate(Math.toRadians(dest.getX()),
Math.toRadians(dest.getY()));
// Get point perpendicular to c2 (c2_90)
gc.setStartingGeographicPoint(Math.toDegrees(c2.x),
Math.toDegrees(c2.y));
gc.setDirection(adjustAngle(c23Dir + 90.0), c23Dist);
dest = gc.getDestinationGeographicPoint();
Coordinate c2_90 = new Coordinate(Math.toRadians(dest.getX()),
Math.toRadians(dest.getY()));
// Find intersecting point from c1-c1_90 and c2-c2_90, this will give us
// one of our pole locations
Coordinate interA = findIntersection(c1, c1_90, c2, c2_90);
// Now intersect great circle of c1-c2 and pole-(aLonR,aLatR). This will
// give us closest/farthest intersecting point
interA = findIntersection(c1, c2, interA, new Coordinate(aLonR, aLatR));
Coordinate interB = new Coordinate(Math.toRadians(adjustAngle(Math
.toDegrees(interA.x + PI))), -interA.y);
// Next we are going to get the angle and distance of interA and interB
// from point c1 to make sure we use the closest point
double radius = getRadius(c1.y);
double[] azm_sidb = azm_sidb(c1.y, c1.x, interA.y, interA.x);
double azimuthA = azm_sidb[0];
double distA = (azm_sidb[1] * radius) / resolution;
azm_sidb = azm_sidb(c1.y, c1.x, interB.y, interB.x);
double azimuthB = azm_sidb[0];
double distB = (azm_sidb[1] * radius) / resolution;
// Now we use the closest intersection point to c1
double azimuth, dist;
Coordinate inter;
if (distA < distB) {
// Closest point is A
azimuth = azimuthA;
dist = distA;
inter = interA;
} else {
// Closest point is B
azimuth = azimuthB;
dist = distB;
inter = interB;
}
// Get angle and distance from intersection point used to initial
// point aLonR, aLatR for bestX value
azm_sidb = azm_sidb(inter.y, inter.x, aLatR, aLonR);
// Convert side_b into projection space
double bestX = (azm_sidb[1] * getRadius(inter.y)) / resolution;
// Figure if bestX is negative or positive
double actual = Math.toRadians(c12Dir);
double expected = azm_sidb[0];
// Correct for world wrapping
if (Math.abs(actual - expected) > Math.PI) {
if (actual > expected) {
expected += TWOPI;
} else {
expected -= TWOPI;
}
}
if (actual < expected) {
bestX = -bestX;
}
// Return {y distance in projection space, difference between angle of
// closest intersecting point and c1 and c1-c2, and our x value in
// projection space}
return new double[] { dist,
Math.abs(adjustAngle(Math.toDegrees(azimuth)) - c12Dir), bestX };
}
/**
* Find the first intersecting point using great circle algorithm for the
* two great circles defined by c1-c2 and c3-c4. Other intersecting point
* will be exactly 180 degrees and -lat on the other side of the world
*
* @param c1
* @param c2
* @param c3
* @param c4
* @return
*/
private static Coordinate findIntersection(Coordinate c1, Coordinate c2,
Coordinate c3, Coordinate c4) {
// Find intersecting points using great circle algorithm...
Coordinate e1Xe2 = correctiveCrossProduct(c1, c2);
Coordinate e3Xe4 = correctiveCrossProduct(c3, c4);
Coordinate ea = normalize(e1Xe2);
Coordinate eb = normalize(e3Xe4);
Coordinate eaXeb = crossProduct(ea, eb);
return new Coordinate(atan2(-eaXeb.y, eaXeb.x), atan2(eaXeb.z,
sqrt(eaXeb.x * eaXeb.x + eaXeb.y * eaXeb.y)));
}
/**
* Adjusts an angle to ensure it is between -180/180
*
* @param angle
* @return
*/
private static double adjustAngle(double angle) {
double newVal = angle % 360;
if (newVal > 180) {
newVal -= 360;
} else if (newVal < -180) {
newVal += 360;
}
return newVal;
}
private static Coordinate correctiveCrossProduct(Coordinate c1,
Coordinate c2) {
Coordinate e1Xe2 = new Coordinate();
double lat1 = c1.y;
double lat2 = c2.y;
double lon1 = c1.x;
double lon2 = c2.x;
e1Xe2.x = sin(lat1 - lat2) * sin((lon1 + lon2) / 2)
* cos((lon1 - lon2) / 2) - sin(lat1 + lat2)
* cos((lon1 + lon2) / 2) * sin((lon1 - lon2) / 2);
e1Xe2.y = sin(lat1 - lat2) * cos((lon1 + lon2) / 2)
* cos((lon1 - lon2) / 2) + sin(lat1 + lat2)
* sin((lon1 + lon2) / 2) * sin((lon1 - lon2) / 2);
e1Xe2.z = cos(lat1) * cos(lat2) * sin(lon1 - lon2);
return e1Xe2;
}
private static Coordinate crossProduct(Coordinate e1, Coordinate e2) {
Coordinate e1Xe2 = new Coordinate();
e1Xe2.x = e1.y * e2.z - e2.y * e1.z;
e1Xe2.y = e1.z * e2.x - e2.z * e1.x;
e1Xe2.z = e1.x * e2.y - e1.y * e2.x;
return e1Xe2;
}
/**
* @param e1Xe2
* @param e1
* @param e2
* @return
*/
private static Coordinate normalize(Coordinate e1Xe2) {
Coordinate ea = new Coordinate();
double length = Math.sqrt(e1Xe2.x * e1Xe2.x + e1Xe2.y * e1Xe2.y
+ e1Xe2.z * e1Xe2.z);
ea.x = e1Xe2.x / length;
ea.y = e1Xe2.y / length;
ea.z = e1Xe2.z / length;
return ea;
}
/*
* (non-Javadoc)
*
@ -820,8 +655,13 @@ public class VIIRSMapProjection extends MapProjection {
}
} else {
side_a = PIO2 - clat;
val_trig = (cos(side_a) * cos(side_b))
+ (sin(side_a) * sin(side_b) * cos(angle_C));
double cos_side_a = cos(side_a);
double sin_side_a = sin(side_a);
double cos_side_b = cos(side_b);
double sin_side_b = sin(side_b);
val_trig = (cos_side_a * cos_side_b)
+ (sin_side_a * sin_side_b * cos(angle_C));
if (val_trig > 1.e0) {
side_c = 0.e0;
} else if (val_trig < -1.e0) {
@ -837,8 +677,8 @@ public class VIIRSMapProjection extends MapProjection {
if ((side_c < 1.e-6) || (side_c == PI)) {
angle_B = 0.e0;
} else {
val_trig = (cos(side_b) - (cos(side_c) * cos(side_a)))
/ (sin(side_c) * sin(side_a));
val_trig = (cos_side_b - (cos(side_c) * cos_side_a))
/ (sin(side_c) * sin_side_a);
if (val_trig > 1.e0) {
angle_B = 0.e0;
} else if (val_trig < -1.e0) {
@ -859,7 +699,7 @@ public class VIIRSMapProjection extends MapProjection {
}
}
return new double[] { Math.toDegrees(alon), Math.toDegrees(alat) };
return new double[] { alon, alat };
}
public static double getRadius(double latitude) {
@ -920,10 +760,16 @@ public class VIIRSMapProjection extends MapProjection {
if (angle_B < -PI)
angle_B = angle_B + TWOPI;
side_b = acos((cos(side_c) * cos(side_a))
+ (sin(side_c) * sin(side_a) * cos(angle_B)));
cos_angle_C = (cos(side_c) - (cos(side_a) * cos(side_b)))
/ (sin(side_a) * sin(side_b));
double cos_side_c = cos(side_c);
double cos_side_a = cos(side_a);
double sin_side_c = sin(side_c);
double sin_side_a = sin(side_a);
double acos_in = (cos_side_c * cos_side_a)
+ (sin_side_c * sin_side_a * cos(angle_B));
side_b = acos(acos_in);
cos_angle_C = (cos_side_c - (cos_side_a * cos(side_b)))
/ (sin_side_a * sin(side_b));
if (cos_angle_C < -9.9999999e-1) {
angle_C = PI;
} else if (cos_angle_C > 9.9999999e-1) {
@ -943,6 +789,14 @@ public class VIIRSMapProjection extends MapProjection {
return new double[] { azimuth, side_b };
}
private static double acos(double a) {
double atan = atan(sqrt(1 - a * a) / a);
if (atan < 0) {
atan += PI;
}
return atan;
}
public static class Provider extends AbstractProvider {
private static final long serialVersionUID = 1266944341235566642L;

View file

@ -980,8 +980,7 @@ public class NctextDbQuery {
}
if(pdc==null){
try {
pdc = PointDataRequest.requestPointDataAllLevels(
null, tableName,
pdc = PointDataRequest.requestPointDataAllLevels(tableName,
parameters.toArray(new String[parameters.size()]),
null, rcMap);
} catch (VizException e) {

View file

@ -544,9 +544,9 @@ public class NcPlotModelDataRequestJob extends Job {
// Datacube didn't have proper plugin; going directly
// to the data store
tempPdc = PointDataRequest.requestPointDataAllLevels(
null, this.plugin,
params.toArray( new String[params.size()] ),
null, map );
this.plugin,
params.toArray(new String[params.size()]), null,
map);
}
long t002 = System.currentTimeMillis();

View file

@ -568,7 +568,7 @@ public class PlotModelGenerator2 extends Job {
// Datacube didn't have proper plugin; going directly
// to the data store
tempPdc = PointDataRequest.requestPointDataAllLevels(
null, this.plugin,
this.plugin,
params.toArray(new String[params.size()]),
null, map);
}