CADRE-old/example_gs.py at master · thearn/CADRE-old · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
from openmdao.lib.casehandlers.api import CSVCaseRecorder
from openmdao.main.api import Component
from openmdao.main.datatypes.api import Float
from CADRE.CADRE_mdp import CADRE_Optimization
import urllib
import json
import numpy as np


class Elevation(Component):

    """
    Computes elevation with respect to mean sea level on Earth at a given
    latitude and longitude.

    Uses google maps API.
    """

    lon = Float(-83.7264, iotype="in", units="deg", descr="Longitude")
    lat = Float(42.2708, iotype="in", units="deg", descr="Latitude")
    alt = Float(0.256, iotype="out", units="km", descr="Elevation")

    lon_ = Float(-83.7264, iotype="out")
    lat_ = Float(42.2708, iotype="out")

    def __init__(self):
        super(Elevation, self).__init__()
        self.J = np.zeros(2)
        self.resolution = 0.001

    def execute(self):
        self.lat_ = self.lat
        self.lon_ = self.lon
        lt1 = self.lat
        ln1 = self.lon
        lt2 = self.lat + self.resolution
        ln2 = self.lon
        lt3 = self.lat
        ln3 = self.lon + self.resolution
        try:
            url = "http://maps.googleapis.com/maps/api/elevation/json?locations=%s,%s|%s,%s|%s,%s&sensor=false" % (
                str(lt1), str(ln1), str(lt2), str(ln2), str(lt3), str(ln3))
            response = ''.join(urllib.urlopen(url).readlines())
            data = json.loads(response)

            self.alt = float(data["results"][0]["elevation"]) / 1000.

            e1 = float(data["results"][1]["elevation"]) / 1000.
            e2 = float(data["results"][2]["elevation"]) / 1000.

            self.J[0] = (e1 - self.alt) / self.resolution
            self.J[1] = (e2 - self.alt) / self.resolution
        except:
            pass

    def linearize(self):
        pass

    def apply_deriv(self, arg, result):
        if "alt" in result:
            if "lat" in arg:
                result["arg"] += self.J[0] * arg["lat"]
            if "lon" in arg:
                result["arg"] += self.J[1] * arg["lon"]

    def apply_derivT(self, arg, result):
        if "alt" in arg:
            if "lat" in result:
                result["lat"] += self.J[0] * arg["alt"]
            if "lon" in result:
                result["lon"] += self.J[1] * arg["alt"]

top = CADRE_Optimization(n=1500, m=300)
top.add("Elevation", Elevation())
top.driver.workflow.add("Elevation")
for i in xrange(6):
    top.connect("Elevation.alt", "pt%s.alt" % str(i))
lons = ["pt" + str(i) + ".lon" for i in xrange(6)] + ["Elevation.lon"]
lats = ["pt" + str(i) + ".lat" for i in xrange(6)] + ["Elevation.lat"]
top.driver.add_parameter(lons, low=-180, high=180)
top.driver.add_parameter(lats, low=-90, high=90)

top.driver.recorders = [CSVCaseRecorder(filename='CADRE_gs.csv')]
printvars = []
for var in ['Data', 'ConCh', 'ConDs', 'ConS0', 'ConS1', 'SOC']:
    printvars += ["pt" + str(i) + ".Data" for i in xrange(6)]
top.driver.printvars = printvars + ["Elevation.alt"]
top.run()