debian-jplephem/jplephem/daf.py

259 lines
9.9 KiB
Python

"""Access a NASA JPL SPICE Double Precision Array File (DAF).
http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html
"""
import io
import mmap
import sys
from struct import Struct
from numpy import array as numpy_array, ndarray
FTPSTR = b'FTPSTR:\r:\n:\r\n:\r\x00:\x81:\x10\xce:ENDFTP' # FTP test string
LOCFMT = {b'BIG-IEEE': '>', b'LTL-IEEE': '<'}
K = 1024
class DAF(object):
"""Access to NASA SPICE Double Precision Array Files (DAF).
Provide the constructor with a ``file_object`` for full access to
both the segment summaries and to the numeric arrays. If you pass a
``StringIO`` instead, then you can fetch the summary information but
not access the arrays.
"""
def __init__(self, file_object):
if getattr(file_object, 'encoding', None):
raise ValueError('file_object must be opened in binary "b" mode')
self.file = file_object
self._map = None
self._array = None
file_record = self.read_record(1)
def unpack():
fmt = self.endian + '8sII60sIII8s603s28s297s'
self.file_record_struct = Struct(fmt)
(locidw, self.nd, self.ni, self.locifn, self.fward, self.bward,
self.free, locfmt, self.prenul, self.ftpstr, self.pstnul
) = self.file_record_struct.unpack(file_record)
self.locidw = file_record[:8].upper().rstrip()
if self.locidw == b'NAIF/DAF':
for self.locfmt, self.endian in LOCFMT.items():
unpack()
if self.nd == 2:
break
else:
raise ValueError('neither a big- nor a little-endian scan'
' of this file produces the expected ND=2')
elif self.locidw.startswith(b'DAF/'):
if file_record[500:1000].strip(b'\0') != FTPSTR:
raise ValueError('this SPK file has been damaged')
self.locfmt = file_record[88:96]
self.endian = LOCFMT.get(self.locfmt)
if self.endian is None:
raise ValueError('unknown format {0!r}'.format(self.locfmt))
unpack()
else:
raise ValueError('file starts with {0!r}, not "NAIF/DAF" or "DAF/"'
.format(self.locidw))
self.locifn_text = self.locifn.rstrip()
summary_format = 'd' * self.nd + 'i' * self.ni
self.summary_control_struct = Struct(self.endian + 'ddd')
self.summary_struct = struct = Struct(self.endian + summary_format)
self.summary_length = length = struct.size
self.summary_step = length + (-length % 8) # pad to 8 bytes
self.summaries_per_record = (1024 - 8 * 3) // self.summary_step
def read_record(self, n):
"""Return record `n` as 1,024 bytes; records are indexed from 1."""
self.file.seek(n * K - K)
return self.file.read(K)
def write_record(self, n, data):
"""Write `data` to file record `n`; records are indexed from 1."""
self.file.seek(n * K - K)
return self.file.write(data)
def write_file_record(self):
data = self.file_record_struct.pack(
self.locidw.ljust(8, b' '), self.nd, self.ni, self.locifn,
self.fward, self.bward, self.free, self.locfmt,
self.prenul, self.ftpstr, self.pstnul,
)
self.write_record(1, data)
def map_words(self, start, end):
"""Return a memory-map of the elements `start` through `end`.
The memory map will offer the 8-byte double-precision floats
("elements") in the file from index `start` through to the index
`end`, inclusive, both counting the first float as element 1.
Memory maps must begin on a page boundary, so `skip` returns the
number of extra bytes at the beginning of the return value.
"""
i, j = 8 * start - 8, 8 * end
try:
fileno = self.file.fileno()
except (AttributeError, io.UnsupportedOperation):
fileno = None
if fileno is None:
skip = 0
self.file.seek(i)
m = self.file.read(j - i)
else:
skip = i % mmap.ALLOCATIONGRANULARITY
r = mmap.ACCESS_READ
m = mmap.mmap(fileno, length=j-i+skip, access=r, offset=i-skip)
if sys.version_info > (3,):
m = memoryview(m) # so further slicing can return views
return m, skip
def comments(self):
"""Return the text inside the comment area of the file."""
record_numbers = range(2, self.fward)
if not record_numbers:
return ''
data = b''.join(self.read_record(n)[0:1000] for n in record_numbers)
try:
return data[:data.find(b'\4')].decode('ascii').replace('\0', '\n')
except IndexError:
raise ValueError('DAF file comment area is missing its EOT byte')
except UnicodeDecodeError:
raise ValueError('DAF file comment area is not ASCII text')
def read_array(self, start, end):
"""Return floats from `start` to `end` inclusive, indexed from 1.
The entire range of floats is immediately read into memory from
the file, making this efficient for small sequences of floats
whose values are all needed immediately.
"""
f = self.file
f.seek(8 * (start - 1))
length = 1 + end - start
data = f.read(8 * length)
return ndarray(length, self.endian + 'd', data)
def map_array(self, start, end):
"""Return floats from `start` to `end` inclusive, indexed from 1.
Instead of pausing to load all of the floats into RAM, this
routine creates a memory map which will load data from the file
only as it is accessed, and then will let it expire back out to
disk later. This is very efficient for large data sets to which
you need random access.
"""
if self._array is None:
self._map, skip = self.map_words(1, self.free - 1)
assert skip == 0
self._array = ndarray(self.free - 1, self.endian + 'd', self._map)
return self._array[start - 1 : end]
def summary_records(self):
"""Yield (record_number, n_summaries, record_data) for each record.
Readers will only use the second two values in each tuple.
Writers can update the record using the `record_number`.
"""
record_number = self.fward
unpack = self.summary_control_struct.unpack
while record_number:
data = self.read_record(record_number)
next_number, previous_number, n_summaries = unpack(data[:24])
yield record_number, n_summaries, data
record_number = int(next_number)
def summaries(self):
"""Yield (name, (value, value, ...)) for each summary in the file."""
length = self.summary_length
step = self.summary_step
for record_number, n_summaries, summary_data in self.summary_records():
name_data = self.read_record(record_number + 1)
for i in range(0, int(n_summaries) * step, step):
j = self.summary_control_struct.size + i
name = name_data[i:i+step].strip()
data = summary_data[j:j+length]
values = self.summary_struct.unpack(data)
yield name, values
def map(self, summary_values):
"""Return the array of floats described by a summary.
Instead of pausing to load all of the floats into RAM, this
routine creates a memory map which will load data from the file
only as it is accessed, and then will let it expire back out to
disk later. This is very efficient for large data sets to which
you need random access.
"""
return self.map_array(summary_values[-2], summary_values[-1])
def add_array(self, name, values, array):
"""Add a new array to the DAF file.
The summary will be initialized with the `name` and `values`,
and will have its start word and end word fields set to point to
where the `array` of floats has been appended to the file.
"""
f = self.file
scs = self.summary_control_struct
record_number = self.bward
data = bytearray(self.read_record(record_number))
next_record, previous_record, n_summaries = scs.unpack(data[:24])
if n_summaries < self.summaries_per_record:
summary_record = record_number
name_record = summary_record + 1
data[:24] = scs.pack(next_record, previous_record, n_summaries + 1)
self.write_record(summary_record, data)
else:
summary_record = ((self.free - 1) * 8 + 1023) // 1024 + 1
name_record = summary_record + 1
free_record = summary_record + 2
n_summaries = 0
data[:24] = scs.pack(summary_record, previous_record, n_summaries)
self.write_record(record_number, data)
summaries = scs.pack(0, record_number, 1).ljust(1024, b'\0')
names = b'\0' * 1024
self.write_record(summary_record, summaries)
self.write_record(name_record, names)
self.bward = summary_record
self.free = (free_record - 1) * 1024 // 8 + 1
start_word = self.free
f.seek((start_word - 1) * 8)
array = numpy_array(array) # TODO: force correct endian
f.write(array.view())
end_word = f.tell() // 8
self.free = end_word + 1
self.write_file_record()
values = values[:self.nd + self.ni - 2] + (start_word, end_word)
base = 1024 * (summary_record - 1)
offset = int(n_summaries) * self.summary_step
f.seek(base + scs.size + offset)
f.write(self.summary_struct.pack(*values))
f.seek(base + 1024 + offset)
f.write(name[:self.summary_length].ljust(self.summary_step, b' '))
NAIF_DAF = DAF # a separate class supported NAIF/DAF format in jplephem 2.2