File: soda.py

package info (click to toggle)
gavodachs 2.3%2Bdfsg-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 7,260 kB
  • sloc: python: 58,359; xml: 8,882; javascript: 3,453; ansic: 661; sh: 158; makefile: 22
file content (281 lines) | stat: -rw-r--r-- 8,972 bytes parent folder | download
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
"""
Helper functions for SODA manipulators.

This primarily comprises various WCS helpers.  It is built on base.coords,
which is where you'll get the wcsFields.

Note that this must not use things from protocols.datalink, as it
is imported from there.  Essentially, use this space for helpers
for SODA manipulations that are generic enough to be kept outside of
the RD but not generic enough to go do base.coords.
"""

#c Copyright 2008-2020, the GAVO project
#c
#c This program is free software, covered by the GNU GPL.  See the
#c COPYING file in the source distribution.


from gavo import base
from gavo import rscdef
from gavo import svcs
from gavo.base import coords
from gavo.utils import fitstools

from gavo.base import makeStruct as MS


DEFAULT_SEMANTICS = "http://dc.g-vo.org/datalink#other"


class EmptyData(base.ExecutiveAction):
	"""raise this when you notice you won't have any data to return.
	"""
	responseCode = 204
	# SODA and HTTP say no bytes are allowed in empty responses.
	responsePayload = ""


class DatalinkFault(object):
	"""A datalink error ("fault", as it's called in the spec).

	These are usually constructed using one of the classmethods

	* AuthenticationFault -- Not authenticated (and authentication required)
	* AuthorizationFault -- Not authorized (to access the resource)
	* NotFoundFault -- Unknown ID value
	* UsageFault -- Invalid input (e.g. no ID values)
	* TransientFault -- Service is not currently able to function
	* FatalFault -- Service cannot perform requested action
	* Fault -- General error (not covered above)

	all of which take the pubDID that caused the failure and a human-oriented
	error message.
	"""
	def __init__(self, code, pubDID, message, exceptionClass, semantics,
			description=None):
		self.code, self.pubDID, self.message = code, pubDID, message
		self.semantics = semantics
		self.exceptionClass = exceptionClass
		self.description = description
	
	@classmethod
	def _addErrorMaker(cls, errCode, exceptionClass):
		def meth(inner, pubDID, message, semantics=DEFAULT_SEMANTICS,
				description=None):
			return inner(errCode, pubDID, message, exceptionClass, semantics,
				description)
		setattr(cls, errCode, classmethod(meth))

	def asDict(self):
		"""returns an error row for the datalink response.
		"""
		return {"ID": self.pubDID, "error_message":
			"%s: %s"%(self.code, self.message),
			"semantics": self.semantics,
			"description": self.description}

	def raiseException(self):
		raise self.exceptionClass(self.message+" (pubDID: %s)"%self.pubDID)

for errName, exClass in [
		("AuthenticationFault", svcs.ForbiddenURI), 
		("AuthorizationFault", svcs.ForbiddenURI),
		("NotFoundFault", svcs.UnknownURI),
		("UsageFault", svcs.BadMethod),
		("TransientFault", svcs.BadMethod),
		("FatalFault", svcs.Error),
		("Fault", svcs.Error)]:
	DatalinkFault._addErrorMaker(errName, exClass)
del errName, exClass


class FormatNow(base.ExecutiveAction):
	"""can be raised by data functions to abort all further processing
	and format the current descriptor.data.
	"""


class DeliverNow(base.ExecutiveAction):
	"""can be raised by data functions to abort all further processing
	and return the current descriptor.data to the client.
	"""


def ensureSkyWCS(descriptor):
	"""furnishes a soda fits descriptor with skyWCS and spatialAxes attributes 
	if necessary.

	This is usually called by the functions needing this automatically, and
	it's a no-op if it has already run.

	skyWCS will be None after this function as run if no usable WCS is found;
	otherwise, it will be a astropy.wcs.WCS instance.  Additionally, spatialAxes
	will be a sequence of 1-based axis indices, and an empty dictionary
	axisNames is available to be filled by metaMakers.  It will then map
	the SODA parameter name to either a FITS axis index or to the special
	WCSLONG, WCSLAT values.
	"""
	if hasattr(descriptor, "skyWCS"):
		return
	
	descriptor.skyWCS, descriptor.spatialAxes = coords.getSkyWCS(descriptor.hdr)
	descriptor.axisNames = {}


def iterSpatialAxisKeys(descriptor, axisMetaOverrides):
	"""yields SODA inputKeys for spatial cutouts along the spatial
	coordinate axes.

	This can be nothing if descriptor doesn't have a skyWCS attribute
	or if it's None.
	"""
	ensureSkyWCS(descriptor)
	if descriptor.skyWCS is None:
		return

	footprint = descriptor.skyWCS.calcFootprint(descriptor.hdr)
	wcsprm = descriptor.skyWCS.wcs

	# FIXME: UCD inference!
	for name, colInd, description, baseUCD, cutoutName in [
		(wcsprm.lattyp.strip(), wcsprm.lat, "The latitude coordinate",
			"pos.eq.dec", "WCSLAT"),
		(wcsprm.lngtyp.strip(), wcsprm.lng, "The longitude coordinate",
			"pos.eq.ra", "WCSLONG")]:
		if name:
			vertexCoos = footprint[:,colInd]
			paramArgs = {"name": name, "unit": "deg", 
					"description": description,
					"ucd": baseUCD}

			minCoo, maxCoo = min(vertexCoos), max(vertexCoos)
			# for RA, we need to move the stitching line out
			# of the way (and go to negative longitudes) if
			# 0 is on the image; we're doing a little heuristic
			# there assuming that images are smaller than 180 deg.
			if cutoutName=="WCSLONG":
				if coords.straddlesStitchingLine(minCoo, maxCoo):
					minCoo, maxCoo = maxCoo-360, minCoo

			if name in axisMetaOverrides:
				paramArgs.update(axisMetaOverrides[name])

			yield MS(svcs.InputKey,  multiplicity="single",
				type="double precision[2]", xtype="interval",
				values=MS(rscdef.Values, min=minCoo, max=maxCoo),
				**paramArgs)
			descriptor.axisNames[name] = cutoutName


def iterOtherAxisKeys(descriptor, axisMetaOverrides):
	"""yields inputKeys for all non-spatial WCS axes.

	descriptor must be a FITSDescriptor.
	"""
	ensureSkyWCS(descriptor)
	if descriptor.skyWCS is None:
		return

	axesLengths = fitstools.getAxisLengths(descriptor.hdr)
	for axIndex, length in enumerate(axesLengths):
		fitsAxis = axIndex+1
		if fitsAxis in descriptor.spatialAxes:
			continue
		if length==1:
			# no cutouts along degenerate axes
			continue
		
		try:
			ax = fitstools.WCSAxis.fromHeader(descriptor.hdr, fitsAxis)
		except ValueError:
			# probably botched WCS, or an inseparable axis.
			# Just ignore this axis, operators can add it manually
			# using forceSeparable
			continue

		descriptor.axisNames[ax.name] = fitsAxis
		minPhys, maxPhys = ax.getLimits()

		# FIXME: ucd inference
		paramArgs = {"name": ax.name, "unit": ax.cunit, 
			"description": "Coordinate along axis number %s"%fitsAxis,
			"ucd": None}
		if fitsAxis in axisMetaOverrides:
			paramArgs.update(axisMetaOverrides[fitsAxis])

		yield MS(svcs.InputKey,  multiplicity="single",
			type="double precision[2]", xtype="interval",
			values=MS(rscdef.Values, min=minPhys, max=maxPhys),
			**paramArgs)


def addPolygonSlices(descriptor, poly, srcPar="Unknown"):
	"""adds slicings in descriptor.slices for a pgsphere.SPoly poly.

	srcPar is the name of the parameter that generated the polygon
	(for making error messages)
	"""
	for axisInd, lower, upper in coords.getPixelLimits(
			poly.asCooPairs(), descriptor.skyWCS):
		descriptor.changingAxis(axisInd, srcPar)
		descriptor.slices.append((axisInd, lower, upper))


def doAxisCutout(descriptor, args):
	"""updates descriptor.data on a FITS descriptor, interpreting the
	parameters defined by iter*AxisKeys, passed in in args.

	This is the main implementation of //soda#fits_doWCSCutout
	"""
	ensureSkyWCS(descriptor)
	slices = descriptor.slices

	# limits: [minRA, maxRA], [minDec, maxDec]]
	footprint = descriptor.skyWCS.calcFootprint(descriptor.hdr)
	limits = [[min(footprint[:,0]), max(footprint[:,0])],
		[min(footprint[:,1]), max(footprint[:,1])]]
	if coords.straddlesStitchingLine(limits[0][0], limits[0][1]):
		limits[0] = [limits[0][1]-360, limits[0][0]]
	limitsChangedName = None

	for parName, fitsAxis in descriptor.axisNames.items():
		if args[parName] is None:
			continue
		limitsChangedName = parName

		if not isinstance(fitsAxis, int):
			# some sort of spherical axis
			if fitsAxis=="WCSLAT":
				cooLimits = limits[1]
			elif fitsAxis=="WCSLONG":
				cooLimits = limits[0]
			else:
				assert False

			cooLimits[0] = max(cooLimits[0], args[parName][0])
			cooLimits[1] = min(cooLimits[1], args[parName][1])
			
		else:
			# 1-d axis
			transform = fitstools.WCSAxis.fromHeader(descriptor.hdr, fitsAxis)
			axMin, axMax = args[parName]
			descriptor.changingAxis(fitsAxis, parName)
			slices.append((fitsAxis, 
				transform.physToPix(axMin), transform.physToPix(axMax)))

	if limitsChangedName:
		for axisInd, lower, upper in coords.getPixelLimits([
				(limits[0][0], limits[1][0]),
				(limits[0][1], limits[1][1])], descriptor.skyWCS):
			descriptor.changingAxis(axisInd, limitsChangedName)
			slices.append((axisInd, lower, upper))

	if slices:
		for axis, lower, upper in slices:
			if lower==upper:  # Sentinel for emtpy data
				raise EmptyData()
		descriptor.data[descriptor.imageExtind] = fitstools.cutoutFITS(
			descriptor.data[descriptor.imageExtind],
			*slices)
		descriptor.dataIsPristine = False