From b9f255274487afbbe44260f955c6283fbd8e9b4a Mon Sep 17 00:00:00 2001 From: AWIPS User Date: Tue, 19 May 2015 17:05:46 -0500 Subject: [PATCH] GVAR McIDAS native projection support added for Raytheon satellite.mcidas decoder. Former-commit-id: bf979e7040709e25b3e8d5262587286a9713aefb --- .../META-INF/MANIFEST.MF | 3 + .../util/satellite/SatSpatialFactory.java | 90 ++++++++++++- .../dataplugin/satellite/SatMapCoverage.java | 30 +++-- .../META-INF/MANIFEST.MF | 3 +- .../mcidas/McidasSatelliteDecoder.java | 126 ++++++++++++++---- 5 files changed, 210 insertions(+), 42 deletions(-) diff --git a/edexOsgi/com.raytheon.edex.plugin.satellite/META-INF/MANIFEST.MF b/edexOsgi/com.raytheon.edex.plugin.satellite/META-INF/MANIFEST.MF index 9eb9e5ea29..6329481bb6 100644 --- a/edexOsgi/com.raytheon.edex.plugin.satellite/META-INF/MANIFEST.MF +++ b/edexOsgi/com.raytheon.edex.plugin.satellite/META-INF/MANIFEST.MF @@ -23,3 +23,6 @@ Require-Bundle: com.raytheon.uf.common.dataplugin;bundle-version="1.12.1174", com.raytheon.uf.edex.database;bundle-version="1.0.0", com.raytheon.uf.edex.menus;bundle-version="1.0.0", com.raytheon.uf.common.numeric;bundle-version="1.14.0" +Import-Package: gov.noaa.nws.ncep.common.tools, + gov.noaa.nws.ncep.edex.util, + org.apache.commons.codec.binary diff --git a/edexOsgi/com.raytheon.edex.plugin.satellite/src/com/raytheon/edex/util/satellite/SatSpatialFactory.java b/edexOsgi/com.raytheon.edex.plugin.satellite/src/com/raytheon/edex/util/satellite/SatSpatialFactory.java index a55f065719..f3b098201c 100644 --- a/edexOsgi/com.raytheon.edex.plugin.satellite/src/com/raytheon/edex/util/satellite/SatSpatialFactory.java +++ b/edexOsgi/com.raytheon.edex.plugin.satellite/src/com/raytheon/edex/util/satellite/SatSpatialFactory.java @@ -20,15 +20,19 @@ package com.raytheon.edex.util.satellite; +import java.awt.geom.Rectangle2D; +import gov.noaa.nws.ncep.common.tools.IDecoderConstantsN; + import org.geotools.geometry.DirectPosition2D; import org.opengis.referencing.crs.ProjectedCRS; import org.opengis.referencing.operation.MathTransform; - import com.raytheon.edex.plugin.satellite.SatelliteDecoderException; import com.raytheon.edex.plugin.satellite.dao.SatMapCoverageDao; import com.raytheon.uf.common.dataplugin.satellite.SatMapCoverage; import com.raytheon.uf.common.geospatial.MapUtil; import com.vividsolutions.jts.geom.Envelope; +import com.vividsolutions.jts.geom.Geometry; +import com.vividsolutions.jts.io.WKTReader; /** * @@ -48,7 +52,7 @@ import com.vividsolutions.jts.geom.Envelope; * Nov 05, 2014 2714 bclement replaced DecoderException with SatelliteDecoderException * Nov 05, 2014 3788 bsteffen use getOrCreateCoverage in place of queryByMapId * May 11, 2015 mjames South polar stereogrpahic support added. - * + * 05/19/2015 mjames@ucar Added decoding of GVAR native projection products * * */ @@ -67,6 +71,8 @@ public class SatSpatialFactory { public static final int PROJ_CYLIN_EQUIDISTANT = 7; public static final int UNDEFINED = -1; + + public static final int PROJ_GVAR = 7585; /** The singleton instance */ private static SatSpatialFactory instance; @@ -210,6 +216,58 @@ public class SatSpatialFactory { return getCoverageSingleCorner(crsType, nx, ny, lov, latin, latin, la1, lo1, dx, dy); } + + public SatMapCoverage getCoverageNative(int crsType, int nx, int ny, + double reflon, int upperLeftElement, int upperLeftLine, + int xres, int yres, ProjectedCRS crs) + throws SatelliteDecoderException { + try { + + // Construct the polygon constructor String + StringBuffer buffer = new StringBuffer(); + buffer.append("POLYGON(("); + buffer.append(reflon - 90. + " -90.0,"); + buffer.append(reflon + 90. + " -90.0,"); + buffer.append(reflon + 90. + " 90.0,"); + buffer.append(reflon - 90. + " 90.0,"); + buffer.append(reflon - 90. + " -90.0"); + buffer.append("))"); + // Create the geometry from the constructed String + Geometry geometry = new WKTReader().read(buffer.toString()); + + // Construct rectangle + double minX = upperLeftElement; + int maxX = upperLeftElement + (nx * xres); + double minY = upperLeftLine + (ny * yres); + minY = -minY; + int maxY = -1 * upperLeftLine; + Rectangle2D rect = new Rectangle2D.Double(minX, + minY, maxX, maxY); + + SatMapCoverage coverage = createCoverageFromNative(crsType, nx, ny, + reflon, upperLeftElement, upperLeftLine, + xres, yres, crs, geometry ); + + return checkPersisted(coverage); + + } catch (Exception e) { + StringBuilder buf = new StringBuilder(); + buf.append( + "Error getting or constructing SatMapCoverage for values: ") + .append("\n\t"); + buf.append("crsType=" + crsType).append("\n\t"); + buf.append("nx=" + nx).append("\n\t"); + buf.append("ny=" + ny).append("\n\t"); + buf.append("reflon=" + reflon).append("\n\t"); + buf.append("upperLeftElement=" + upperLeftElement).append("\n\t"); + buf.append("upperLeftLine=" + upperLeftLine).append("\n\t"); + buf.append("xres=" + xres).append("\n\t"); + buf.append("yres=" + yres).append("\n\t"); + throw new SatelliteDecoderException(buf.toString(), e); + } + } + + /** * * Create a {@link SatMapCoverage} with an area defined by two corners. The @@ -309,6 +367,29 @@ public class SatSpatialFactory { return new SatMapCoverage(crsType, envelope.getMinX(), envelope.getMinY(), nx, ny, dx, dy, crs); } + + /** + * Create a SatMapCoverage from native projection + */ + private static SatMapCoverage createCoverageFromNative(Integer crsType, + Integer nx, Integer ny, double reflon, int upperLeftElement, + int upperLeftLine, int xres, int yres, ProjectedCRS crs, + Geometry geometry) { + + // do your shit here. + float dx = IDecoderConstantsN.FLOAT_MISSING; + float dy = IDecoderConstantsN.FLOAT_MISSING; + + double minX, minY; + //if (crsType == PROJ_GVAR) { // for native projection + minX = upperLeftElement; + minY = upperLeftLine + (ny * yres); + minY = -minY; + + return new SatMapCoverage(crsType, minX, minY, nx, ny, + dx, dy, crs, geometry); + + } /** * Create a {@link ProjectedCRS} from a crsType and some parameters. @@ -345,7 +426,7 @@ public class SatSpatialFactory { case PROJ_LAMBERT: return createLambertCrs(latin, latin2, lov); case PROJ_CYLIN_EQUIDISTANT: - return createEqCylCrs(latin, lov); + return createEqCylCrs(latin, lov); default: return createNorthPolarStereoCrs(latin, lov); } @@ -376,5 +457,8 @@ public class SatSpatialFactory { return MapUtil.constructSouthPolarStereo(MapUtil.AWIPS_EARTH_RADIUS, MapUtil.AWIPS_EARTH_RADIUS, latin, lov); } + + + } diff --git a/edexOsgi/com.raytheon.uf.common.dataplugin.satellite/src/com/raytheon/uf/common/dataplugin/satellite/SatMapCoverage.java b/edexOsgi/com.raytheon.uf.common.dataplugin.satellite/src/com/raytheon/uf/common/dataplugin/satellite/SatMapCoverage.java index 80c0339777..f42c3fe1ce 100644 --- a/edexOsgi/com.raytheon.uf.common.dataplugin.satellite/src/com/raytheon/uf/common/dataplugin/satellite/SatMapCoverage.java +++ b/edexOsgi/com.raytheon.uf.common.dataplugin.satellite/src/com/raytheon/uf/common/dataplugin/satellite/SatMapCoverage.java @@ -80,7 +80,8 @@ import com.vividsolutions.jts.geom.Polygon; * Apr 11, 2014 2947 bsteffen Fix equals * Oct 16, 2014 3454 bphillip Upgrading to Hibernate 4 * Nov 05, 2014 3788 bsteffen Make gid a sequence instead of a hash. - * + * May 19, 2015 mjames@ucar Added decoding of GVAR native projection products, + * increased crsWKT to 5120 for GVAR the_geom * */ @Entity @@ -92,6 +93,8 @@ import com.vividsolutions.jts.geom.Polygon; public class SatMapCoverage extends PersistableDataObject implements IGridGeometryProvider { + public static final int PROJ_GVAR = 7585; + private static final long serialVersionUID = 1L; @Id @@ -102,7 +105,8 @@ public class SatMapCoverage extends PersistableDataObject implements /** * The projection of the map coverage 1=Mercator, 3=Lambert Conformal - * 5=Polar Stereographic + * 5=Polar Stereographic, 7585 = native satellite navigation e.g. + * GVAR, ... * * @deprecated This field is only useful for GINI satellite format decoding * and should not be in the coverage object @@ -110,7 +114,6 @@ public class SatMapCoverage extends PersistableDataObject implements @Column @XmlAttribute @DynamicSerializeElement - @Deprecated private Integer projection; /** Minimum x coordinate in crs space */ @@ -149,11 +152,11 @@ public class SatMapCoverage extends PersistableDataObject implements @DynamicSerializeElement private double dy; - @Column(length = 2047) + @Column(length = 5120) @XmlAttribute @DynamicSerializeElement private String crsWKT; - + @Transient private CoordinateReferenceSystem crsObject; @@ -316,7 +319,7 @@ public class SatMapCoverage extends PersistableDataObject implements public void setDy(double dy) { this.dy = dy; } - + public String getCrsWKT() { if (crsWKT == null && crsObject != null) { crsWKT = crsObject.toWKT(); @@ -379,12 +382,19 @@ public class SatMapCoverage extends PersistableDataObject implements @Override public GridGeometry2D getGridGeometry() { - GridEnvelope gridRange = new GridEnvelope2D(0, 0, getNx(), getNy()); - Envelope crsRange = new Envelope2D(getCrs(), new Rectangle2D.Double( - minX, minY, getNx() * getDx(), getNy() * getDy())); + GridEnvelope gridRange; + Envelope crsRange; + + int nx = getNx(); + int ny = getNy(); + gridRange = new GridEnvelope2D(0, 0, getNx(), getNy()); + crsRange = new Envelope2D(getCrs(), new Rectangle2D.Double( + minX, minY, getNx() * getDx(), getNy() * getDy())); + + return new GridGeometry2D(gridRange, crsRange); } - + @Override public int hashCode() { HashCodeBuilder builder = new HashCodeBuilder(); diff --git a/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/META-INF/MANIFEST.MF b/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/META-INF/MANIFEST.MF index bd61ea5495..d6433c7784 100644 --- a/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/META-INF/MANIFEST.MF +++ b/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/META-INF/MANIFEST.MF @@ -16,5 +16,6 @@ Bundle-RequiredExecutionEnvironment: JavaSE-1.7 Import-Package: com.raytheon.uf.common.localization, com.raytheon.uf.common.menus, com.raytheon.uf.common.menus.xml, - com.raytheon.uf.common.status + com.raytheon.uf.common.status, + org.apache.commons.codec.binary Export-Package: com.raytheon.uf.edex.plugin.satellite.mcidas diff --git a/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/src/com/raytheon/uf/edex/plugin/satellite/mcidas/McidasSatelliteDecoder.java b/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/src/com/raytheon/uf/edex/plugin/satellite/mcidas/McidasSatelliteDecoder.java index 09f0bf1114..848f88a6af 100644 --- a/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/src/com/raytheon/uf/edex/plugin/satellite/mcidas/McidasSatelliteDecoder.java +++ b/edexOsgi/com.raytheon.uf.edex.plugin.satellite.mcidas/src/com/raytheon/uf/edex/plugin/satellite/mcidas/McidasSatelliteDecoder.java @@ -26,6 +26,9 @@ import java.util.Calendar; import java.util.GregorianCalendar; import java.util.TimeZone; +import org.apache.commons.codec.binary.Base64; +import org.opengis.referencing.crs.ProjectedCRS; + import com.raytheon.edex.esb.Headers; import com.raytheon.edex.exception.DecoderException; import com.raytheon.edex.util.satellite.SatSpatialFactory; @@ -33,6 +36,7 @@ import com.raytheon.uf.common.dataplugin.PluginDataObject; import com.raytheon.uf.common.dataplugin.satellite.SatMapCoverage; import com.raytheon.uf.common.dataplugin.satellite.SatelliteRecord; import com.raytheon.uf.common.datastorage.records.IDataRecord; +import com.raytheon.uf.common.geospatial.MapUtil; import com.raytheon.uf.common.status.IUFStatusHandler; import com.raytheon.uf.common.status.UFStatus; import com.raytheon.uf.common.time.DataTime; @@ -71,7 +75,8 @@ import com.raytheon.uf.edex.plugin.satellite.mcidas.util.McidasSatelliteLookups. * IDataRecord required by the SatelliteDao * 12/03/2013 DR 16841 D. Friedman Allow record overwrites * 09/18/2014 3627 mapeters Updated deprecated method calls. - * 05/11/2015 mjames PS (south and north) stereogrpahic support added. + * 05/11/2015 mjames@ucar PS (south and north) stereogrpahic support added. + * 05/19/2015 mjames@ucar Added decoding of GVAR native projection products * * * @author @@ -90,6 +95,8 @@ public class McidasSatelliteDecoder { private static final int RADIUS = 6371200; + final int SIZE_OF_AREA = 256; + private static final double HALFPI = Math.PI / 2.; private static final double RTD = 180. / Math.PI; @@ -134,9 +141,17 @@ public class McidasSatelliteDecoder { * @throws Exception */ private PluginDataObject[] decodeMcidasArea(byte[] data) throws Exception { + + + byte[] area = null; + byte[] nonAreaBlock = new byte[data.length - SIZE_OF_AREA]; + area = new byte[SIZE_OF_AREA]; + System.arraycopy(data, 0, area, 0, SIZE_OF_AREA); + System.arraycopy(data, SIZE_OF_AREA, nonAreaBlock, 0, + nonAreaBlock.length); + ByteBuffer buf = ByteBuffer.wrap(data); buf.order(ByteOrder.LITTLE_ENDIAN); - // Decode the directory block if (buf.getInt() != 0) { formatError(UNEXPECTED_HEADER_VALUE); @@ -148,24 +163,37 @@ public class McidasSatelliteDecoder { formatError(UNEXPECTED_HEADER_VALUE); } } - int sensorSourceNumber = buf.getInt(); - int yyyddd = buf.getInt(); - int hhmmss = buf.getInt(); - int ulImageLine = buf.getInt(); - int ulImageElement = buf.getInt(); - buf.getInt(); // reserved - int nLines = buf.getInt(); - int nElementsPerLine = buf.getInt(); - int nBytesPerElement = buf.getInt(); - int lineResolution = buf.getInt(); - int elementResolution = buf.getInt(); - int nBands = buf.getInt(); - int linePrefixLength = buf.getInt(); - /* int projectNumber = */buf.getInt(); - /* int creationYyyddd = */buf.getInt(); - /* int creationHhmmss = */buf.getInt(); - int bandMap1to32 = buf.getInt(); - int bandMap33to64 = buf.getInt(); + + int sensorSourceNumber = buf.getInt(); // W3 + int yyyddd = buf.getInt(); // W4 + int hhmmss = buf.getInt(); // W5 + int ulImageLine = buf.getInt(); // W6 + int ulImageElement = buf.getInt(); // W7 + buf.getInt(); // reserved // W8 + int nLines = buf.getInt(); // W9 + int nElementsPerLine = buf.getInt(); // W10 + int nBytesPerElement = buf.getInt(); // W11 + int lineResolution = buf.getInt(); // W12 + int elementResolution = buf.getInt(); // W13 + int nBands = buf.getInt(); // W14 + int linePrefixLength = buf.getInt(); // W15 + /* int projectNumber = */buf.getInt(); // W16 + /* int creationYyyddd = */buf.getInt(); // W17 + /* int creationHhmmss = */buf.getInt(); // W18 + + /* + * W19 + 32-bit filter band map for multichannel + images; if a bit is set, data exists for the band; + band 1 is the least significant byte (rightmost) + */ + int bandMap1to32 = buf.getInt(); // W19 + /* + * W20-24 + satellite specific information + */ + int bandMap33to64 = buf.getInt(); // W20 + buf.position(buf.position() + (4 * 4)); // sensor specific buf.position(buf.position() + (4 * 8)); // memo int areaNumber = buf.getInt(); @@ -188,12 +216,35 @@ public class McidasSatelliteDecoder { /* int scaling = */buf.getInt(); /* int supplementalBlockOffset = */buf.getInt(); buf.getInt(); // reserved - /* int calibrationOffset = */buf.getInt(); + int calibrationOffset = buf.getInt(); buf.getInt(); // comment cards - + + + int navsize; + if (calibrationOffset == 0){ + navsize = dataBlockOffset - navBlockOffset; + } else { + navsize = calibrationOffset - navBlockOffset; + } + byte[] navigation = new byte[navsize]; + System.arraycopy(nonAreaBlock, 0, navigation, 0, navsize); + + /* mjames@ucar + * bandMap1to32 is a 32-bit filter band map for multichannel images. + * if a bit is set, data exists for the band; band 1 is the least + * significant byte (rightmost) + * + * Example: for GVAR GEWCOMP UNIWISC image, + * bandMap1to32 = 4 + * bandMap33to64 = -1, so + * nBands = 1 + * bandBitsCount = 33 + * + * this is a problem... + */ long bandBits = ((long) bandMap33to64 << 32) | bandMap1to32; long bandBitsCount = Long.bitCount(bandBits); - if (nBands != bandBitsCount) { + if (nBands != bandBitsCount && nBands > 1) { formatError("Specified number of bands does not match number of bits in band map"); } @@ -201,7 +252,7 @@ public class McidasSatelliteDecoder { buf.position(navBlockOffset); SatMapCoverage coverage = decodeNavigation(elementResolution, lineResolution, ulImageElement, ulImageLine, nElementsPerLine, - nLines, buf); + nLines, buf, navigation); // Decode the data block, creating a SatelliteRecord for each band. PluginDataObject[] result = new PluginDataObject[nBands]; @@ -280,7 +331,8 @@ public class McidasSatelliteDecoder { * */ private SatMapCoverage decodeNavigation(int xImgRes, int yImgRes, int ulX, - int ulY, int nx, int ny, ByteBuffer buf) throws Exception { + int ulY, int nx, int ny, ByteBuffer buf, byte[] navigation) + throws Exception { SatMapCoverage result = new SatMapCoverage(); String navType = get4cc(buf); int lineOfEquator = buf.getInt(); @@ -288,7 +340,8 @@ public class McidasSatelliteDecoder { int stdLatDDMMSS = buf.getInt(); int spacingAtStdLatInMeters = buf.getInt(); int nrmlLonDDMMSS = buf.getInt(); - + + // NOTE: We do not check the following for compatibility with WGS84. int radiusInMeters = buf.getInt(); /* int eccentricity = */buf.getInt(); @@ -373,9 +426,20 @@ public class McidasSatelliteDecoder { result = SatSpatialFactory.getInstance().getCoverageTwoCorners( SatSpatialFactory.PROJ_POLAR, nx, ny, (float) clon, (float) clat, la1, lo1, la2, lo2); - + } else { - unimplemented(String.format("navigation type \"%s\"", navType)); + /* + * Native projection + */ + //unimplemented(String.format("navigation type \"%s\"", navType)); + clon = nrmlLonDDMMSS / 10000000.f; + clon = (float) Math.toDegrees(clon); + ProjectedCRS crs = MapUtil.constructNative(navType, encodeNavBlock(navigation)); + + result = SatSpatialFactory.getInstance().getCoverageNative( + SatSpatialFactory.PROJ_GVAR, nx, ny, (float) clon, + ulX, ulY, xImgRes, yImgRes, crs); + } return result; @@ -451,5 +515,11 @@ public class McidasSatelliteDecoder { throw new DecoderException(String.format("%s: unimplemented: %s", traceId, feature)); } + private String encodeNavBlock(byte[] navigation) { + Base64 b64 = new Base64(); + byte[] coded = b64.encode(navigation); + + return new String(coded); + } }