我正在使用基于NDB的Google App Engine应用程序,该应用程序需要跟踪大量(〜2000个)固定位置的昼/夜周期。由于纬度和经度永远不会改变,因此我可以使用PyEphem之类的东西预先对其进行预先计算。我正在使用NDB。如我所见,可能的策略是:


要将一年中的日出值预先计算为日期时间对象,请将
将它们放入列表中,腌制列表,然后将其放入PickleProperty
,但将列表放入JsonProperty
使用DateTimeProperty并设置repeat = True


现在,我希望索引下一个日出/日落属性,但是可以从列表中弹出该属性并将其放入其自己的DateTimeProperty中,以便我可以定期使用查询来确定哪些位置已更改为其他部分周期的。整个列表不需要索引。

有谁知道相对的努力-就这三种方法的索引编制和CPU负载而言?重复=真对索引有影响吗?

谢谢,
戴夫

最佳答案

答案“仅在实例启动时对其进行计算”或“将这些结构预计算并输出到硬编码的python结构中”的答案似乎忽略了乘以一年的日出值所需要的365倍,或者如果是2000倍则存储实例启动时完成计算。使用pyEphem,计算2000日出和日落需要超过两秒钟的时间。在源代码中为2000个位置存储一年的日出和日落可能会占用20兆字节以上的空间。如果有效地腌制这些数字,则需要2 * 365 * 2000 * 8 = 11,680,000字节。

一种更快更好的方法是为一个位置的时间建立最小二乘模型,以建立另一个位置的时间。如下所述,这使总使用空间减少了约70倍。

首先,如果点A和B在相同的纬度上,并且具有相似的高度和地平线参数,则A处的日出发生在相对于B处日出的恒定时间偏移处。例如,如果A在B处以西15度,则日出发生在在A点比在B点晚一小时。其次,如果点A,B,C在相同的经度和低纬度处,则可以相当准确地将一个点的日出时间计算为另外两个点的线性组合。在高纬度或更高的精度下,可以使用多个时间曲线的线性组合。第三,可以将3月20日A点(spring equinox)的日出时间用作归一化点,因此所有计算都可以归一化为同一纬度。

下表显示了使用四个时间曲线的线性组合可以产生什么样的精度。对于距赤道最大46°的经度,结果保持在大约半秒内。对于48°至60°,结果保持在5秒以内。在64°时,结果可能会长达2分钟,而在65°时,可能会导致约6分钟。但是这些时间可能足以满足大多数实际目的。请注意,以下所示程序在66°时崩溃,因为它无法处理pyEphem引发的异常;即使66°低于Arctic Circle 66.5622°N,也会发生“ AlwaysUpError:'Sun'仍在2013/6/14 07:20:15的地平线上”。

修改程序很容易,以便它使用所需的任意数量的时间曲线(请参阅程序中的各种lata = ...语句),给出所需的精度,但要以存储更多曲线和更多系数为代价。当然,可以改变模型以使用时间曲线的子集。例如,可以存储10条曲线,并根据最接近任何给定目标纬度的4条纬度进行计算。但是,对于此演示程序,尚无此类改进。

Lat.  0.0:  Error range: -0.000000 to 0.000000 seconds
Lat.  5.0:  Error range: -0.370571 to 0.424092 seconds
Lat. 10.0:  Error range: -0.486193 to 0.557997 seconds
Lat. 15.0:  Error range: -0.414288 to 0.477041 seconds
Lat. 20.0:  Error range: -0.213614 to 0.247057 seconds
Lat. 25.0:  Error range: -0.065826 to 0.056358 seconds
Lat. 30.0:  Error range: -0.382425 to 0.323623 seconds
Lat. 35.0:  Error range: -0.585914 to 0.488351 seconds
Lat. 40.0:  Error range: -0.490303 to 0.400563 seconds
Lat. 45.0:  Error range: -0.164706 to 0.207415 seconds
Lat. 47.0:  Error range: -0.590103 to 0.756647 seconds
Lat. 48.0:  Error range: -0.852844 to 1.102608 seconds
Lat. 50.0:  Error range: -1.478688 to 1.940351 seconds
Lat. 55.0:  Error range: -3.342506 to 4.696076 seconds
Lat. 60.0:  Error range: -0.000002 to 0.000003 seconds
Lat. 61.0:  Error range: -7.012057 to 4.273954 seconds
Lat. 62.0:  Error range: -21.374033 to 12.347188 seconds
Lat. 63.0:  Error range: -51.872753 to 27.853411 seconds
Lat. 64.0:  Error range: -124.000365 to 59.661029 seconds
Lat. 65.0:  Error range: -351.425224 to 139.656187 seconds


使用上述方法,对于2000个位置中的每个位置,您需要存储五个浮点数:3月20日的日出时间,以及四个时间曲线的四个乘数系数。 (前面提到的减少70倍是因为每个位置存储5个数字,而不是365个数字。)对于每个时间曲线,将存储365个数字,输入项i是与3月20日的日出时差。存储四个时间曲线所用的空间是存储2000条时间曲线的1/500,因此曲线存储空间以乘数为主导。

在给出使用scipy.optimize.leastsq求解系数的程序之前,这里有两个代码片段,可以在ipython解释器中使用,以制作精度表并绘制图表以可视化错误。

import sunrise as sr
for lat in range(0, 65, 5):
    sr.lsr(lat, -110, 2013, 4)


上面产生了前面显示的大多数错误表。 lsr的第三个参数称为daySkip,值4使lsr每隔四天(即一年中仅约90天)工作一次,以进行更快的测试。使用sr.lsr(lat, -110, 2013, 1)会产生相似的结果,但所需时间是原来的四倍。

sr.plotData(15,1./(24*3600))


上面告诉告诉Sunrise.plotData绘制所有内容(要近似的日出数据;模型的近似结果;残差(按秒缩放)以及基数曲线。)

该程序如下所示。请注意,它主要针对北半球经度进行了测试。如果时间曲线足够对称,则程序将按原样处理南半球的经度;如果误差太大,则可以将南半球经度添加到基本曲线中,或者可以更改模型以在赤道以南使用单独的一组曲线。请注意,日落不是在此程序中计算的。对于日落,添加类似于next_setting(ephem.Sun())调用的previous_rising(ephem.Sun())调用,并存储另外四个时间曲线。

#!/usr/bin/python
import ephem, numpy, scipy, scipy.optimize

# Make a set of observers (observation points)
def observers(lata, lona):
    def makeIter(x):
        if hasattr(x, '__iter__'):
            return x
        return [x]
    lata, lona = makeIter(lata), makeIter(lona)
    arr = []
    for lat in lata:
        for lon in lona:
            o = ephem.Observer()
            o.lat, o.lon, o.elevation, o.temp = str(lat), str(lon), 1400, 0
            arr.append(o)
    return tuple(arr)

# Make a year of data for an observer, equinox-relative
def riseData(observr, year, skip):
    yy = ephem.Date('{} 00:00'.format(year))
    rr = numpy.arange(0.0, 366.0, skip)
    springEquinox = 78
    observr.date = ephem.Date(yy + springEquinox)
    seDelta = observr.previous_rising(ephem.Sun()) - yy - springEquinox + 1
    for i, day in enumerate(range(0, 366, skip)):
        observr.date = ephem.Date(yy + day)
        t = observr.previous_rising(ephem.Sun()) - yy - day + 1 - seDelta
        rr[i] = t
    return numpy.array(rr)

# Make a set of time curves
def makeRarSet(lata, lona, year, daySkip):
    rar=[]
    for o in observers(lata, lona):
        r = riseData(o, year, daySkip)
        rar.append(r)
    x = numpy.arange(0., 366., daySkip)
    return (x, rar)

# data() is an object that stores curves + results
def data(s):
    return data.ss[s]

# Initialize data.ss
def setData(lata, lona, year, daySkip):
    x, rar = makeRarSet(lata, lona, year, daySkip)
    data.ss = rar

# Compute y values from model, given vector x and given params in p
def yModel(x, p):
    t = numpy.zeros(len(x))
    for i in range(len(p)):
        t += p[i] * data(i)
    return t

# Compute residuals, given params in p and data in x, y vectors.
# x = independent var,  y = dependent = observations
def residuals(p, y, x):
    err = y - yModel(x, p)
    return err

# Compute least squares result
def lsr(lat, lon, year, daySkip):
    latStep = 13.
    lata = numpy.arange(0., 66.4, latStep)
    lata = [ 88 * (1 - 1.2**-i) for i in range(8)]
    l, lata, lstep, ldown = 0, [], 20, 3
    l, lata, lstep, ldown = 0, [], 24, 4
    while l < 65:
        lata.append(l); l += lstep; lstep -= ldown
    #print 'lata =', lata
    setData(lata, lon, year, daySkip)
    x, ya = makeRarSet(lat, lon, year, daySkip)
    x, za = makeRarSet(lat, 10+lon, year, daySkip)
    data.ss.append(za[0])
    y = ya[0]
    pini = [(0 if abs(lat-l)>latStep else 0.5) for l in lata]
    pars = scipy.optimize.leastsq(residuals, pini, args=(y, x))
    data.x, data.y, data.pv = x, y, yModel(x, pars[0])
    data.par, data.err = pars, residuals(pars[0], y, x)
    #print 'pars[0] = ', pars[0]
    variance = numpy.inner(data.err, data.err)/len(y)
    #print 'variance:', variance
    sec = 1.0/(24*60*60)
    emin, emax = min(data.err), max(data.err)
    print ('Lat. {:4.1f}:  Error range: {:.6f} to {:.6f} seconds'.format(lat, emin/sec, emax/sec))

def plotData(iopt, emul):
    import matplotlib.pyplot as plt
    plt.clf()
    x = data.x
    if iopt == 0:
        iopt = 15
        emul = 1
    if iopt & 1:
        plt.plot(x, data.y)
        plt.plot(x, data.y + 0.001)
        plt.plot(x, data.y - 0.001)
    if iopt & 2:
        plt.plot(x, data.pv)
    if iopt & 4:
        plt.plot(x, emul*data.err)
    if iopt & 8:
        for ya in data.ss:
            plt.plot(x, ya)
    plt.show()

10-06 08:51