Java tutorial
/* * This file is part of JGrasstools (http://www.jgrasstools.org) * (C) HydroloGIS - www.hydrologis.com * * JGrasstools is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ package org.jgrasstools.hortonmachine.modules.hydrogeomorphology.insolation; import static org.jgrasstools.gears.libs.modules.ModelsEngine.calcInverseSunVector; import static org.jgrasstools.gears.libs.modules.ModelsEngine.calcNormalSunVector; import static org.jgrasstools.gears.libs.modules.ModelsEngine.calculateFactor; import static org.jgrasstools.gears.libs.modules.ModelsEngine.scalarProduct; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_AUTHORCONTACTS; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_AUTHORNAMES; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_KEYWORDS; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_LABEL; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_LICENSE; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_NAME; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_STATUS; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_inElev_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_outIns_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_tEndDate_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSINSOLATION_tStartDate_DESCRIPTION; import java.awt.image.RenderedImage; import java.awt.image.SampleModel; import java.awt.image.WritableRaster; import java.util.HashMap; import javax.media.jai.RasterFactory; import javax.media.jai.iterator.RandomIter; import javax.media.jai.iterator.RandomIterFactory; import javax.media.jai.iterator.WritableRandomIter; import oms3.annotations.Author; import oms3.annotations.Description; import oms3.annotations.Execute; import oms3.annotations.In; import oms3.annotations.Keywords; import oms3.annotations.Label; import oms3.annotations.License; import oms3.annotations.Name; import oms3.annotations.Out; import oms3.annotations.Status; import org.geotools.coverage.grid.GridCoverage2D; import org.geotools.referencing.crs.DefaultGeographicCRS; import org.jgrasstools.gears.libs.modules.JGTModel; import org.jgrasstools.gears.utils.CrsUtilities; import org.jgrasstools.gears.utils.coverage.CoverageUtilities; import org.jgrasstools.gears.utils.geometry.GeometryUtilities; import org.jgrasstools.hortonmachine.i18n.HortonMessageHandler; import org.joda.time.DateTime; import org.joda.time.DateTimeZone; import org.joda.time.format.DateTimeFormat; import org.joda.time.format.DateTimeFormatter; import org.opengis.referencing.crs.CoordinateReferenceSystem; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Point; @Description(OMSINSOLATION_DESCRIPTION) @Author(name = OMSINSOLATION_AUTHORNAMES, contact = OMSINSOLATION_AUTHORCONTACTS) @Keywords(OMSINSOLATION_KEYWORDS) @Label(OMSINSOLATION_LABEL) @Name(OMSINSOLATION_NAME) @Status(OMSINSOLATION_STATUS) @License(OMSINSOLATION_LICENSE) public class OmsInsolation extends JGTModel { @Description(OMSINSOLATION_inElev_DESCRIPTION) @In public GridCoverage2D inElev = null; @Description(OMSINSOLATION_tStartDate_DESCRIPTION) @In public String tStartDate = null; @Description(OMSINSOLATION_tEndDate_DESCRIPTION) @In public String tEndDate = null; @Description(OMSINSOLATION_outIns_DESCRIPTION) @Out public GridCoverage2D outIns; private static final double pCmO3 = 0.3; private static final double pRH = 0.4; private static final double pLapse = -.0065; private static final double pVisibility = 60; /** * The solar constant. */ private static final double SOLARCTE = 1368.0; /** * The atmosphere pressure. */ private static final double ATM = 1013.25; private double lambda; private double delta; private double omega; private HortonMessageHandler msg = HortonMessageHandler.getInstance(); @Execute public void process() throws Exception { // transform the checkNull(inElev, tStartDate, tEndDate); // extract some attributes of the map HashMap<String, Double> attribute = CoverageUtilities.getRegionParamsFromGridCoverage(inElev); double dx = attribute.get(CoverageUtilities.XRES); /* * The models use only one value of the latitude. So I have decided to * set it to the center of the raster. Extract the CRS of the * GridCoverage and transform the value of a WGS84 latitude. */ CoordinateReferenceSystem sourceCRS = inElev.getCoordinateReferenceSystem2D(); CoordinateReferenceSystem targetCRS = DefaultGeographicCRS.WGS84; double srcPts[] = new double[] { attribute.get(CoverageUtilities.EAST), attribute.get(CoverageUtilities.SOUTH) }; Coordinate source = new Coordinate(srcPts[0], srcPts[1]); Point[] so = new Point[] { GeometryUtilities.gf().createPoint(source) }; CrsUtilities.reproject(sourceCRS, targetCRS, so); // the latitude value lambda = Math.toRadians(so[0].getY()); /* * transform the start and end date in an int value (the day in the * year, from 1 to 365) */ DateTimeFormatter formatter = DateTimeFormat.forPattern("yyyy-MM-dd").withZone(DateTimeZone.UTC); DateTime currentDatetime = formatter.parseDateTime(tStartDate); int startDay = currentDatetime.getDayOfYear(); currentDatetime = formatter.parseDateTime(tEndDate); int endDay = currentDatetime.getDayOfYear(); CoverageUtilities.getRegionParamsFromGridCoverage(inElev); RenderedImage pitTmpRI = inElev.getRenderedImage(); int width = pitTmpRI.getWidth(); int height = pitTmpRI.getHeight(); WritableRaster pitWR = CoverageUtilities.replaceNovalue(pitTmpRI, -9999.0); pitTmpRI = null; WritableRaster insolationWR = CoverageUtilities.createDoubleWritableRaster(width, height, null, pitWR.getSampleModel(), 0.0); WritableRandomIter insolationIterator = RandomIterFactory.createWritable(insolationWR, null); WritableRaster gradientWR = normalVector(pitWR, dx); pm.beginTask(msg.message("insolation.calculating"), endDay - startDay); for (int i = startDay; i <= endDay; i++) { calcInsolation(lambda, pitWR, gradientWR, insolationWR, i, dx); pm.worked(i - startDay); } pm.done(); for (int y = 2; y < height - 2; y++) { for (int x = 2; x < width - 2; x++) { if (pitWR.getSampleDouble(x, y, 0) == -9999.0) { insolationIterator.setSample(x, y, 0, Double.NaN); } } } outIns = CoverageUtilities.buildCoverage("insolation", insolationWR, attribute, inElev.getCoordinateReferenceSystem()); } /** * Evaluate the radiation. * * @param lambda * the latitude. * @param demWR * the raster of elevation * @param gradientWR * the raster of the gradient value of the dem. * @param insolationWR * the wr where to store the result. * @param the * day in the year. * @paradx the resolutiono of the dem. */ private void calcInsolation(double lambda, WritableRaster demWR, WritableRaster gradientWR, WritableRaster insolationWR, int day, double dx) { // calculating the day angle // double dayang = 2 * Math.PI * (day - 1) / 365.0; double dayangb = (360 / 365.25) * (day - 79.436); dayangb = Math.toRadians(dayangb); // Evaluate the declination of the sun. delta = getDeclination(dayangb); // Evaluate the radiation in this day. double ss = Math.acos(-Math.tan(delta) * Math.tan(lambda)); double hour = -ss + (Math.PI / 48.0); while (hour <= ss - (Math.PI / 48)) { omega = hour; // calculating the vector related to the sun double sunVector[] = calcSunVector(); double zenith = calcZenith(sunVector[2]); double[] inverseSunVector = calcInverseSunVector(sunVector); double[] normalSunVector = calcNormalSunVector(sunVector); int height = demWR.getHeight(); int width = demWR.getWidth(); WritableRaster sOmbraWR = calculateFactor(height, width, sunVector, inverseSunVector, normalSunVector, demWR, dx); double mr = 1 / (sunVector[2] + 0.15 * Math.pow((93.885 - zenith), (-1.253))); for (int j = 0; j < height; j++) { for (int i = 0; i < width; i++) { // evaluate the radiation. calcRadiation(i, j, demWR, sOmbraWR, insolationWR, sunVector, gradientWR, mr); } } hour = hour + Math.PI / 24.0; } } /* * Evaluate the declination. */ private double getDeclination(double dayangb) { double delta = .3723 + 23.2567 * Math.sin(dayangb) - .758 * Math.cos(dayangb) + .1149 * Math.sin(2 * dayangb) + .3656 * Math.cos(2 * dayangb) - .1712 * Math.sin(3 * dayangb) + .0201 * Math.cos(3 * dayangb); return Math.toRadians(delta); } /* * evaluate several component of the radiation and then multiply by the * sOmbra factor. */ private void calcRadiation(int i, int j, WritableRaster demWR, WritableRaster sOmbraWR, WritableRaster insolationWR, double[] sunVector, WritableRaster gradientWR, double mr) { double z = demWR.getSampleDouble(i, j, 0); double pressure = ATM * Math.exp(-0.0001184 * z); double ma = mr * pressure / ATM; double temp = 273 + pLapse * (z - 4000); double vap_psat = Math.exp(26.23 - 5416.0 / temp); double wPrec = 0.493 * pRH * vap_psat / temp; double taur = Math.exp((-.09030 * Math.pow(ma, 0.84)) * (1.0 + ma - Math.pow(ma, 1.01))); double d = pCmO3 * mr; double tauo = 1 - (0.1611 * d * Math.pow(1.0 + 139.48 * d, -0.3035) - 0.002715 * d) / (1.0 + 0.044 * d + 0.0003 * Math.pow(d, 2)); double taug = Math.exp(-0.0127 * Math.pow(ma, 0.26)); double tauw = 1 - 2.4959 * (wPrec * mr) / (1.0 + 79.034 * (wPrec * mr) * 0.6828 + 6.385 * (wPrec * mr)); double taua = Math.pow((0.97 - 1.265 * Math.pow(pVisibility, (-0.66))), Math.pow(ma, 0.9)); double In = 0.9751 * SOLARCTE * taur * tauo * taug * tauw * taua; double cosinc = scalarProduct(sunVector, gradientWR.getPixel(i, j, new double[3])); if (cosinc < 0) { cosinc = 0; } double tmp = insolationWR.getSampleDouble(i, j, 0); insolationWR.setSample(i, j, 0, In * cosinc * sOmbraWR.getSampleDouble(i, j, 0) / 1000 + tmp); } protected double[] calcSunVector() { double sunVector[] = new double[3]; sunVector[0] = -Math.sin(omega) * Math.cos(delta); sunVector[1] = Math.sin(lambda) * Math.cos(omega) * Math.cos(delta) - Math.cos(lambda) * Math.sin(delta); sunVector[2] = Math.cos(lambda) * Math.cos(omega) * Math.cos(delta) + Math.sin(lambda) * Math.sin(delta); return sunVector; } protected WritableRaster normalVector(WritableRaster pitWR, double res) { int minX = pitWR.getMinX(); int minY = pitWR.getMinY(); int rows = pitWR.getHeight(); int cols = pitWR.getWidth(); RandomIter pitIter = RandomIterFactory.create(pitWR, null); /* * Initializa the Image of the normal vector in the central point of the * cells, which have 3 components so the Image have 3 bands.. */ SampleModel sm = RasterFactory.createBandedSampleModel(5, cols, rows, 3); WritableRaster tmpNormalVectorWR = CoverageUtilities.createDoubleWritableRaster(cols, rows, null, sm, 0.0); WritableRandomIter tmpNormalIter = RandomIterFactory.createWritable(tmpNormalVectorWR, null); /* * apply the corripio's formula (is the formula (3) in the article) */ for (int j = minY; j < minX + rows - 1; j++) { for (int i = minX; i < minX + cols - 1; i++) { double zij = pitIter.getSampleDouble(i, j, 0); double zidxj = pitIter.getSampleDouble(i + 1, j, 0); double zijdy = pitIter.getSampleDouble(i, j + 1, 0); double zidxjdy = pitIter.getSampleDouble(i + 1, j + 1, 0); double firstComponent = res * (zij - zidxj + zijdy - zidxjdy); double secondComponent = res * (zij + zidxj - zijdy - zidxjdy); double thirthComponent = 2 * (res * res); double den = Math.sqrt(firstComponent * firstComponent + secondComponent * secondComponent + thirthComponent * thirthComponent); tmpNormalIter.setPixel(i, j, new double[] { firstComponent / den, secondComponent / den, thirthComponent / den }); } } pitIter.done(); return tmpNormalVectorWR; } private double calcZenith(double sunVector2) { return Math.acos(sunVector2); } }