diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms index 8a5d693a83e22922fb6f08fcac195ab51feec70c..d20564f6df0649777126c12084ccaf41b56e4c11 100644 --- a/GAMS/IntExtOpt.gms +++ b/GAMS/IntExtOpt.gms @@ -61,7 +61,7 @@ $gdxin oilcrops 0.7 pulses 0.5 starchyRoots 4.0 - pasture 0.1 / ; + pasture 0.001 / ; SCALAR fertiliserUnitCost / 1.0 / @@ -69,6 +69,7 @@ $gdxin area(crop, location) total area for each crop - Mha fertI(crop, location) fertilizer intensity for each crop - factor between 0 and 1 irrigI(crop, location) irrigation intensity for each crop - factor between 0 and 1 + otherIntensity(crop, location) feedAmount(crop) amount of feed for each crop - Mt netImportAmount(all_types) net imports of crops and meat - Mt yield(crop, location) yield per area for each crop - t per ha @@ -81,7 +82,7 @@ $gdxin pastureDecrease(location) energy total input energy - energy; - POSITIVE VARIABLE area, fertI, irrigI, feedAmount, agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease; + POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount, agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease; EQUATIONS UNIT_ENERGY_EQ(crop, location) energy per area @@ -92,6 +93,7 @@ $gdxin MEAT_DEMAND_CONSTRAINT satisfy demand for meat MAX_FERT_INTENSITY_CONSTRAINT(crop, location) constraint on maximum fertilizer intensity MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity + MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) TOTAL_LAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change NON_FEED_CROP_CONSTRAINT(not_feed_crop) constraint to set non feed crop feed usage to zero MAX_NET_IMPORT_CONSTRAINT(import_crop) constraint on max net imports @@ -106,9 +108,11 @@ $gdxin PASTURE_DECREASE_CONV_CALC(location) ENERGY_EQ total energy objective function; - UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= ( baseCost(crop) * (cropAdj(crop)) + + UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= ( baseCost(crop) + fertiliserUnitCost * fertI(crop, location) + - irrigCost(location) * irrigI(crop, location) ) ; + irrigCost(location) * irrigI(crop, location) + + SQR(otherIntensity(crop, location)) * 1 + ) ; YIELD_EQ(crop, location) .. yield(crop, location) =E= ( yieldNone(crop, location) + @@ -116,7 +120,7 @@ $gdxin (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) + (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) * ((fertI(crop, location) + delta) ** fertParam(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) - ) * cropAdj(crop); + ) * cropAdj(crop) * otherIntensity(crop, location); NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + netImportAmount(crop); @@ -130,6 +134,7 @@ $gdxin MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1; MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1; + MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) .. otherIntensity(crop, location) =L= 1; TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, area(crop_less_pasture, location)) / (1.0 - unhandledCropArea) + area('pasture', location); @@ -148,10 +153,10 @@ $gdxin PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= area('pasture', location) - previousArea('pasture', location); PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(area('pasture', location) - previousArea('pasture', location)); - ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location)) + ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location)) + (sum(location, agriLandExpansion(location)) * landChangeEnergy) - + (sum(location, SQR((cropIncrease(location) + cropDecrease(location)) / suitableLandArea(location)) * (cropIncrease(location) + cropDecrease(location)) * landChangeEnergy * 10)) - + (sum(location, SQR((pastureIncrease(location) + pastureDecrease(location)) / suitableLandArea(location)) * (pastureIncrease(location) + pastureDecrease(location)) * landChangeEnergy * 10)) + + (sum(location, (0.2 * SQR((cropIncrease(location) + cropDecrease(location)) / suitableLandArea(location)) + (cropIncrease(location) + cropDecrease(location)) * 0.1) * (cropIncrease(location) + cropDecrease(location)) * landChangeEnergy)) + + (sum(location, (0.2 * SQR((pastureIncrease(location) + pastureDecrease(location)) / suitableLandArea(location)) + (pastureIncrease(location) + pastureDecrease(location)) * 0.1) * (pastureIncrease(location) + pastureDecrease(location)) * landChangeEnergy)) + sum(import_crop, (netImportAmount(import_crop)) * worldInputPrices(import_crop))) / 1000000; MODEL LAND_USE /ALL/ ; @@ -169,21 +174,21 @@ $gdxin * if fitToPreviousAreas is zero then we just optimise if (fitToPreviousAreas eq 0, - fertI.L(crop, location) = 0.5; - irrigI.L(crop, location) = 0.5; - feedAmount.L(crop) = 0.0; - netImportAmount.L(import_crop) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2; - area.L(crop, location) = previousArea(crop, location) - SOLVE LAND_USE USING NLP MINIMIZING energy; - - - + fertI.L(crop, location) = 0.5; + irrigI.L(crop, location) = 0.5; + otherIntensity.L(crop, location) = 0.1; + netImportAmount.L(import_crop) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2; + area.L(crop, location) = previousArea(crop, location) + SOLVE LAND_USE USING NLP MINIMIZING energy; + + + * if fitToPreviousAreas is not zero then we try to find a cropAdj factors that give the previous areas ie we attempt to calibrate else basePasture = sum(location, previousArea("pasture", location)); baseCropland = sum((location, crop_less_pasture), previousArea(crop_less_pasture, location)); - while((count le 10 and not (abs(basePasture-newPasture) le 0.1 and abs(baseCropland-newCropland) le 0.1)), + while((count le 10 and not (abs(basePasture-newPasture) le 0.01 and abs(baseCropland-newCropland) le 0.1)), fertI.L(crop, location) = 0.5; irrigI.L(crop, location) = 0.5; area.L(crop, location) = previousArea(crop, location); @@ -209,7 +214,7 @@ $gdxin display area.l, count, basePasture, baseCropland, newPasture, newCropland, cropAdj; count = count+1 ; ); - ); + ); display agriLandExpansion.l, previousArea, area.l, net_supply.l, demand, feedAmount.l, yield.l, netImportAmount.l; diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java index ee041bf37a4f4ed7515b8131060371dff88652a6..d22f1d5a37d5b20041e1ebf7ea215e816932196f 100644 --- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java +++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java @@ -39,13 +39,17 @@ public class GamsCountryInput { this.calibrateToObserved = calibrateToObserved; } + public GamsCountryInput(GamsCountryInput gamsInput, Map<CropType, Double> cropAdjs) { + this(gamsInput.country, gamsInput.projectedDemand, gamsInput.worldInputPrices, gamsInput.maxNetImport, gamsInput.minNetImport, cropAdjs, false); + } + public static GamsCountryInput createInput(CompositeCountry country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputPrices, Map<CropType, Double> baseNetImport, Map<CropType, Double> maxOfProdOrSupply, Map<CropType, Double> cropAdjustments, boolean calibrateToObserved) { Map<CropType, Double> maxNetImport; Map<CropType, Double> minNetImport; - if (calibrateToObserved) { + if (calibrateToObserved) { // this condition is no required as calibration now done from GasmRasterOptimiser rather than in GAMS script maxNetImport = baseNetImport; minNetImport = baseNetImport; diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java index baca8ee9317fb80e87e5a4814f8d3404af18ca2b..f780d39a714d235aa24b7445998c8c13bad46b20 100644 --- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java +++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java @@ -5,6 +5,7 @@ import java.util.Map; import ac.ed.lurg.Timestep; import ac.ed.lurg.landuse.AreasItem; import ac.ed.lurg.landuse.IrrigationCostItem; +import ac.ed.lurg.types.CropType; import ac.ed.lurg.yield.YieldResponsesItem; public class GamsLocationInput { @@ -25,6 +26,11 @@ public class GamsLocationInput { this.countryInput = countryInput; } + public static GamsLocationInput createWithNewAdjustments(GamsLocationInput inputToCopy, Map<CropType, Double> cropAdjustments) { + GamsCountryInput updatedGamsCountryInput = new GamsCountryInput(inputToCopy.getCountryInput(), cropAdjustments); + return new GamsLocationInput(inputToCopy.timestep, inputToCopy.yields, inputToCopy.previousAreas, inputToCopy.irrigationCosts, updatedGamsCountryInput); + } + public Map<Integer, ? extends YieldResponsesItem> getYields() { return yields; } diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java index 7b7d24d8058d3c4a8d3d797f90b6b535c6d11b43..fa044bb5d343a05f3d9ea55fe3e4fe5b20d536e1 100644 --- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java @@ -205,7 +205,7 @@ public class GamsLocationOptimiser { GAMSVariable varAreas = outDB.getVariable("area"); GAMSVariable varFertIntensities = outDB.getVariable("fertI"); GAMSVariable varIrrigIntensities = outDB.getVariable("irrigI"); - // GAMSVariable varOtherIntensities = outDB.getVariable("otherIntensity"); + GAMSVariable varOtherIntensities = outDB.getVariable("otherIntensity"); GAMSVariable varFeedAmount = outDB.getVariable("feedAmount"); GAMSVariable varNetImports = outDB.getVariable("netImportAmount"); GAMSVariable varYields = outDB.getVariable("yield"); @@ -215,7 +215,8 @@ public class GamsLocationOptimiser { GAMSParameter parmProdCost = outDB.getParameter("totalProdCost"); GAMSParameter parmCroplandArea = outDB.getParameter("totalCropland"); - double totalArea = 0; + double totalCropArea = 0; + double totalPastureArea = 0; double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, feedAmount, netImport, yield, unitEnergy, cropAdj, prod, prodCost; LazyHashMap<Integer, IntensitiesItem> intensities = new LazyHashMap<Integer, IntensitiesItem>() { @@ -234,7 +235,7 @@ public class GamsLocationOptimiser { area = rec.getLevel(); fertIntensity = varFertIntensities.findRecord(itemName, locationName).getLevel(); irrigIntensity = varIrrigIntensities.findRecord(itemName, locationName).getLevel(); - // otherIntensity = varOtherIntensities.findRecord(itemName, locationName).getLevel(); + otherIntensity = varOtherIntensities.findRecord(itemName, locationName).getLevel(); yield = varYields.findRecord(itemName, locationName).getLevel(); unitEnergy = varUnitEnergies.findRecord(itemName, locationName).getLevel(); @@ -261,13 +262,16 @@ public class GamsLocationOptimiser { double croplandArea = getParmValue(parmCroplandArea, locationName); AreasItem areasItem = cropAreas.lazyGet(locId); - if (cropType.equals(CropType.PASTURE)) + if (cropType.equals(CropType.PASTURE)) { areasItem.setLandCoverArea(LandCoverType.PASTURE, area); - else + totalPastureArea += area; + } + else { areasItem.setLandCoverArea(LandCoverType.CROPLAND, croplandArea); // will set this multiple times, once for each arable crop, but doesn't really matter - + totalCropArea += area; + } + areasItem.setCropFraction(cropType, croplandArea > 0 ? area/croplandArea : 0); - totalArea += area; } netImport = varNetImports.findRecord(CropType.MEAT.getGamsName()).getLevel(); @@ -277,7 +281,7 @@ public class GamsLocationOptimiser { cropUsageData.put(CropType.MEAT, new CropUsageData(0.0, netImport, prod, prodCost)); if (DEBUG) { LogWriter.println(String.format("\n%s:\t\t\t\t\tnetImports= %.3f,\tprod= %.3f,\tprodCost= %.3f", CropType.MEAT.getGamsName(), netImport, prod, prodCost)); - LogWriter.println(String.format("\nTotal area= %.1f", totalArea)); + LogWriter.println(String.format("\nTotal area= %.1f (crop=%.1f, pasture %.1f)", totalCropArea+totalPastureArea, totalCropArea, totalPastureArea)); } GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, cropUsageData, cropAdjs); diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java index a35d7298b40d4c55d7c04e04d89319b9d09ad962..963dcd1460319cb71901932410ab14121088b7fe 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java @@ -2,6 +2,7 @@ package ac.ed.lurg.country.gams; import java.util.ArrayList; import java.util.Collections; +import java.util.HashMap; import java.util.HashSet; import java.util.List; import java.util.Map; @@ -38,14 +39,64 @@ public class GamsRasterOptimiser { // workout similar areas GamsLocationInput gamsInput = convertFromRaster(rasterInputData); - // run optimizer - GamsLocationOptimiser opti = new GamsLocationOptimiser(gamsInput); - GamsLocationOutput gamsOutput = opti.run(); + long start = System.currentTimeMillis(); + + GamsLocationOutput gamsOutput = runGamsWithCalib(gamsInput); + + long last = System.currentTimeMillis(); + LogWriter.printlnError("Took " + (last-start)); // map results back to raster return convertToRaster(gamsInput, gamsOutput); } + /** Probably should be in GAMS class */ + private GamsLocationOutput runGamsWithCalib(GamsLocationInput gamsInput) { + + GamsLocationOutput result = null; + + // this code does similar calibration to that in the GAMS script + if (gamsInput.getCountryInput().isCalibrateToObserved()) { + Map<CropType, Double> cropAdjustments = new HashMap<CropType, Double>(); // recreate, but might have been null is this case anyway + for (CropType c : CropType.getNonMeatTypes()) { + cropAdjustments.put(c, 1.0); + } + + double baseCropland = AreasItem.getTotalLandCover(gamsInput.getPreviousAreas(), LandCoverType.CROPLAND); + double basePasture = AreasItem.getTotalLandCover(gamsInput.getPreviousAreas(), LandCoverType.PASTURE); + + LogWriter.println("baseCropland=" + baseCropland + ", basePasture=" + basePasture); + + for (int i=0; i<10; i++) { + + GamsLocationInput newInput = GamsLocationInput.createWithNewAdjustments(gamsInput, cropAdjustments); + GamsLocationOptimiser opti= new GamsLocationOptimiser(newInput); + result = opti.run(); + + if (baseCropland > 0) { + double cropAdjSingle = AreasItem.getTotalLandCover(result.getCropAreas(), LandCoverType.CROPLAND) / baseCropland; + cropAdjSingle = Math.min(Math.max(cropAdjSingle, 0.5), 2); + + for (CropType c : CropType.getCropsLessPasture()) { + cropAdjustments.put(c, Math.min(Math.max(cropAdjustments.get(c) * cropAdjSingle, 0.2), 5)); + } + } + + if (basePasture > 0) { + double pastureAdjSingle = AreasItem.getTotalLandCover(result.getCropAreas(), LandCoverType.PASTURE) / basePasture; + pastureAdjSingle = Math.min(Math.max(pastureAdjSingle, 0.5), 2); + cropAdjustments.put(CropType.PASTURE, Math.min(Math.max(cropAdjustments.get(CropType.PASTURE) * pastureAdjSingle, 0.2), 5)); + } + } + } + else { + GamsLocationOptimiser opti = new GamsLocationOptimiser(gamsInput); + result = opti.run(); + } + return result; + } + + private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) { RasterSet<AreasItem> newAreaRaster = allocAreas(gamsInput.getPreviousAreas(), gamsOutput); RasterSet<IntensitiesItem> newIntensityRaster = allocIntensities(gamsOutput); @@ -56,33 +107,33 @@ public class GamsRasterOptimiser { for (RasterKey key : entry.getValue()) locationIdRaster.put(key, new IntegerRasterItem(locId)); } - + return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData(), gamsOutput.getCropAdjustments(), locationIdRaster); } private RasterSet<AreasItem> createWithSameLandCovers(RasterSet<AreasItem> toCopy) { RasterSet<AreasItem> theCopy = new RasterSet<AreasItem>(toCopy.getHeaderDetails()); - + for (Entry<RasterKey, AreasItem> entry : toCopy.entrySet()) { AreasItem newAreasItem = new AreasItem(); newAreasItem.setLandCoverAreas(entry.getValue()); theCopy.put(entry.getKey(), newAreasItem); } - + return theCopy; } - + private void checkedTotalAreas(RasterSet<AreasItem> areaRaster, String comment) { for (LandCoverType l : LandCoverType.values()) { double total = 0; for (AreasItem a : areaRaster.values()) { total += a.getLandCoverArea(l); } - + LogWriter.printlnError("Total Area " + comment + ": " + l.getName() + ": " + total); } } - + private RasterSet<AreasItem> allocAreas(Map<Integer, ? extends AreasItem> prevAreasAgg, GamsLocationOutput gamsOutput) { RasterSet<AreasItem> newAreaRaster = createWithSameLandCovers(rasterInputData.getPreviousAreas()); @@ -97,7 +148,7 @@ public class GamsRasterOptimiser { if (DEBUG) checkedTotalAreas(newAreaRaster.popSubsetForKeys(new RasterSet<AreasItem>(newAreaRaster.getHeaderDetails()), keys), locId + " before"); - + double pastureChange = newAreaAggItem.getLandCoverArea(LandCoverType.PASTURE) - prevAreaAggItem.getLandCoverArea(LandCoverType.PASTURE); double croplandChange = newAreaAggItem.getLandCoverArea(LandCoverType.CROPLAND) - prevAreaAggItem.getLandCoverArea(LandCoverType.CROPLAND); double prevForestToNaturalFraction = prevAreaAggItem.getForestToNaturalFraction(); @@ -108,7 +159,7 @@ public class GamsRasterOptimiser { if (pastureChange > 0 && croplandChange < 0) { pastureFromCrop = Math.min(pastureChange, -croplandChange); - + if (pastureChange > -croplandChange) pastureFromNatural = pastureChange+croplandChange; else @@ -116,7 +167,7 @@ public class GamsRasterOptimiser { } else if (pastureChange < 0 && croplandChange > 0) { pastureFromCrop = -Math.min(-pastureChange, croplandChange); - + if (croplandChange > -pastureChange) cropFromNatural = pastureChange+croplandChange; else @@ -126,7 +177,7 @@ public class GamsRasterOptimiser { pastureFromNatural = pastureChange; cropFromNatural = croplandChange; } - + double shortfall = 0; shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.PASTURE, LandCoverType.FOREST, pastureFromNatural * prevForestToNaturalFraction); shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.PASTURE, LandCoverType.OTHER_NATURAL, pastureFromNatural * (1-prevForestToNaturalFraction)); @@ -150,18 +201,18 @@ public class GamsRasterOptimiser { return newAreaRaster; } - + private double allocAllLandCrops(RasterSet<AreasItem> newAreaRaster, Set<RasterKey> keys, LandCoverType toLC, LandCoverType fromLC, double change) { LandCoverType toType = toLC; LandCoverType fromType = fromLC; - + // reverse direction if negative if (change < 0) { change = -change; toType = fromLC; fromType = toLC; } - + double totalShortfall = 0; Set<RasterKey> keysWithSpace = new HashSet<RasterKey>(); double avgChange = change/keys.size(); @@ -169,7 +220,7 @@ public class GamsRasterOptimiser { for (RasterKey key : keys) { //if (DEBUG) LogWriter.println(" Processing raster key " + key); AreasItem newAreasRasterItem = newAreaRaster.get(key); - + double shortfall = newAreasRasterItem.moveAreas(toType, fromType, avgChange); if (shortfall == 0) keysWithSpace.add(key); @@ -205,8 +256,8 @@ public class GamsRasterOptimiser { // as a first attempt only going to look at pasture and cereal yields, assuming other yields will be correlated to one or the other. YieldRaster yieldRaster = rasterInputData.getYields(); RasterSet<AreasItem> cropAreaRaster = rasterInputData.getPreviousAreas(); - - /* checkedTotalAreas(cropAreaRaster, "before"); + + /* checkedTotalAreas(cropAreaRaster, "before"); new RasterOutputer<AreasItem>(cropAreaRaster, "BeforeLandCover") { @Override public Double getValue(RasterKey location) { @@ -237,8 +288,8 @@ public class GamsRasterOptimiser { RasterSet<IrrigationCostItem> irrigCostRaster = rasterInputData.getIrrigationCost(); { - // YieldResponsesItem yresp = yieldRaster.get(new RasterKey(40,38));//getFromCoordinates(-118.0, -33.0); - // LogWriter.printlnError("Test key2: " + yresp.getYieldMax(CropType.WHEAT) + ", " + yresp.getYieldFertOnly(CropType.WHEAT) + ", " + yresp.getYieldIrrigOnly(CropType.WHEAT)); + // YieldResponsesItem yresp = yieldRaster.get(new RasterKey(40,38));//getFromCoordinates(-118.0, -33.0); + // LogWriter.printlnError("Test key2: " + yresp.getYieldMax(CropType.WHEAT) + ", " + yresp.getYieldFertOnly(CropType.WHEAT) + ", " + yresp.getYieldIrrigOnly(CropType.WHEAT)); } // look for inconsistencies @@ -280,8 +331,8 @@ public class GamsRasterOptimiser { int countFound = 0, countMissing = 0; - - + + for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) { YieldResponsesItem yresp = entry.getValue(); RasterKey key = entry.getKey(); @@ -299,7 +350,7 @@ public class GamsRasterOptimiser { LogWriter.printlnError("GamsRasterOptimiser: Can't find cropAreas for " + key); continue; } - + IrrigationCostItem irrigCost = irrigCostRaster.get(key); int cerealCat = findCategory(wheatlDivisions, yresp.getYieldNone(CropType.WHEAT) + yresp.getYieldMax(CropType.WHEAT)); @@ -328,11 +379,11 @@ public class GamsRasterOptimiser { double areaSoFar = aggAreas.getCropArea(crop); double areaThisTime = cropAreas.getCropArea(crop); double updateCropland = (cropAreas.getLandCoverArea(LandCoverType.CROPLAND) + aggAreas.getLandCoverArea(LandCoverType.CROPLAND)); - + if (updateCropland > 0) aggAreas.setCropFraction(crop, (areaSoFar + areaThisTime) / updateCropland); } - + // Land covers ares for (LandCoverType landType : LandCoverType.values()) { double areaThisTime = cropAreas.getLandCoverArea(landType); @@ -346,8 +397,8 @@ public class GamsRasterOptimiser { for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) { LogWriter.println(e.getKey() + " zone has " + e.getValue().size() + " raster areas"); - - CropType[] cs = {CropType.WHEAT, CropType.MAIZE}; + + /* CropType[] cs = {CropType.WHEAT, CropType.MAIZE}; for (CropType c : cs) { for (RasterKey key : e.getValue()) { YieldResponsesItem yri = yieldRaster.get(key); @@ -356,7 +407,7 @@ public class GamsRasterOptimiser { yri.getYieldNone(c), yri.getYieldMidFertOnly(c), yri.getYieldFertOnly(c), yri.getYieldIrrigOnly(c), yri.getYieldMax(c))); } LogWriter.println(""); - } + } */ } return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, rasterInputData.getCountryInput()); diff --git a/src/ac/ed/lurg/landuse/AreasItem.java b/src/ac/ed/lurg/landuse/AreasItem.java index 01b0181deed0542085fb6c0d0cfe0a8b8738b15c..2056c57074fd8a46ca9a8b2b152cf19455b73cef 100644 --- a/src/ac/ed/lurg/landuse/AreasItem.java +++ b/src/ac/ed/lurg/landuse/AreasItem.java @@ -148,4 +148,13 @@ public class AreasItem implements InterpolatingRasterItem<AreasItem> { Double res = from + factor * (to - from); return res; } + + public static double getTotalLandCover(Map<?, ? extends AreasItem> areaRaster, LandCoverType landCover) { + double total = 0; + for (AreasItem a : areaRaster.values()) { + total += a.getLandCoverArea(landCover); + } + + return total; + } }