eoxserver.resources.coverages.domainset

1 #------------------------------------------------------------------------------- 2 # $Id: domainset.py 2240 2013-02-19 19:56:56Z schindlerf $ 3 # 4 # Project: EOxServer <http://eoxserver.org> 5 # Authors: Stephan Krause <stephan.krause@eox.at> 6 # Stephan Meissl <stephan.meissl@eox.at> 7 # 8 #------------------------------------------------------------------------------- 9 # Copyright (C) 2011 EOX IT Services GmbH 10 # 11 # Permission is hereby granted, free of charge, to any person obtaining a copy 12 # of this software and associated documentation files (the "Software"), to deal 13 # in the Software without restriction, including without limitation the rights 14 # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 15 # copies of the Software, and to permit persons to whom the Software is 16 # furnished to do so, subject to the following conditions: 17 # 18 # The above copyright notice and this permission notice shall be included in all 19 # copies of this Software or works derived from this Software. 20 # 21 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 22 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 23 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 24 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 25 # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 26 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 27 # THE SOFTWARE. 28 #------------------------------------------------------------------------------- 29 30 import logging 31 import numpy.linalg 32 33 from django.contrib.gis.geos import GEOSGeometry, Polygon, LineString 34 35 from eoxserver.core.exceptions import InternalError 36 from eoxserver.core.util.timetools import getDateTime 37 from eoxserver.services.exceptions import ( 38 InvalidAxisLabelException, InvalidSubsettingException 39 ) 40 41 42 logger = logging.getLogger(__name__) 43

44 -class RectifiedGrid(object):

45 - def __init__(self, 46 dim=2, 47 spatial_dim=2, 48 low=(0,0), 49 high=(0,0), 50 axis_labels=('lon','lat'), 51 srid=4326, 52 origin=(0,0), 53 offsets=((1,0),(0,1)) 54 ):

55 super(RectifiedGrid, self).__init__() 56 self.dim = dim 57 self.spatial_dim = spatial_dim 58 self.low = low 59 self.high = high 60 self.axis_labels = axis_labels 61 self.srid = srid 62 self.origin = origin 63 self.offsets = offsets

64 65 # TODO: validate inputs 66

67 - def getExtent2D(self):

68 if self.spatial_dim >= 2: 69 lowx = self.origin[0] + self.low[0] * self.offsets[0][0] + self.low[1] * self.offsets[1][0] 70 lowy = self.origin[1] + self.low[0] * self.offsets[0][1] + self.low[1] * self.offsets[1][1] 71 highx = self.origin[0] + self.high[0] * self.offsets[0][0] + self.high[1] * self.offsets[1][0] 72 highy = self.origin[1] + self.high[0] * self.offsets[0][1] + self.high[1] * self.offsets[1][1] 73 74 return (min(lowx, highx), min(lowy, highy), max(lowx, highx), max(lowy, highy)) 75 else: 76 raise Exception("Cannot compute 2D extent of grid with less than 2 spatial dimensions")

77

78 - def getBBOX(self):

79 bbox = Polygon.from_bbox(self.getExtent2D()) 80 bbox.srid = int(self.srid) 81 82 return bbox

83

84 - def contains(self, grid):

85 this_minx, this_miny, this_maxx, this_maxy = self.getExtent2D() 86 that_minx, that_miny, that_maxx, that_maxy = grid.getExtent2D() 87 88 return this_minx <= that_minx and this_miny <= that_miny and \ 89 this_maxx >= that_maxx and this_maxy >= that_maxy

90

91 - def isSubGrid(self, grid):

92 logger.debug("EOxSRectifiedGrid.isSubGrid: Own Extent: (%f,%f,%f,%f)" % self.getExtent2D()) 93 logger.debug("EOxSRectifiedGrid.isSubGrid: Containing Extent: (%f,%f,%f,%f)" % grid.getExtent2D()) 94 95 if grid.contains(self): 96 this_offsets = numpy.array([self.offsets[i][0:2] for i in range(0, 2)]) 97 that_offsets = numpy.array([grid.offsets[i][0:2] for i in range(0, 2)]) 98 99 if not numpy.all(numpy.equal(numpy.abs(this_offsets), numpy.abs(that_offsets))): 100 return False 101 else: 102 v = numpy.linalg.solve(this_offsets, numpy.array(self.origin[0:2]) - numpy.array(grid.origin[0:2])) 103 104 logger.debug("EOxSRectifiedGrid.isSubGrid: v=(%f,%f)" % (v[0], v[1])) 105 106 EPSILON = 1e-10 107 108 if numpy.all(numpy.less(numpy.absolute((numpy.rint(v) - v) / numpy.maximum(v, 1)), EPSILON)): 109 return True 110 else: 111 return False 112 else: 113 return False

114

115 -class Subsetting(object):

116 - def __init__(self, dimension, crs):

117 super(Subsetting, self).__init__() 118 119 self.dimension = dimension 120 self.crs = crs

121

122 - def normalize(self, dimension, value):

123 if value is None or len(value) == 0: 124 return None 125 elif dimension in ("phenomenonTime", "time", "t"): 126 if value[0] == '"' and value[-1] == '"': 127 token = value.lstrip('"').rstrip('"') 128 return getDateTime(token) # this raises an UnknkownParameterFormatException if the datetime format is not recognized 129 else: 130 raise InvalidSubsettingException("Date/Time tokens have to be enclosed in quotation marks (\")") 131 else: 132 try: 133 return float(value) 134 except: 135 raise InvalidSubsettingException("'%s' not recognized as a number" % value)

136

137 - def validate(self, grid=None):

138 return True

139

140 - def _getDataFromFootprint(self, footprint):

141 srid = 4326 142 143 env_minx, env_miny, env_maxx, env_maxy = footprint.extent 144 145 return (srid, env_minx, env_miny, env_maxx, env_maxy)

146

147 -class Slice(Subsetting):

148 - def __init__(self, dimension, crs, slice_point):

149 super(Slice, self).__init__(dimension, crs) 150 151 self.slice_point = self.normalize(dimension, slice_point)

152

153 - def validate(self, grid=None):

154 if self.slice_point is not None: 155 if grid is not None: 156 if self.dimension not in grid.axis_labels: 157 raise InvalidAxisLabelException("Invalid axis label '%s'" % self.dimension) 158 else: 159 raise InvalidSubsettingException("Empty slices are not allowed") 160 161 return True

162

163 - def crosses(self, footprint):

164 srid, env_minx, env_miny, env_maxx, env_maxy = self._getDataFromFootprint(footprint) 165 166 if self.dimension in ("long", "Long"): 167 line = LineString((self.slice_point, env_miny), (self.slice_point, env_maxy), srid=srid) 168 elif self.dimension in ("lat", "Lat"): 169 line = LineString((env_minx, self.slice_point), (env_maxx, self.slice_point), srid=srid) 170 else: 171 raise InternalError("Can handle 2D coverages only.") 172 173 # TODO bbox not defined 174 #return line.crosses(bbox) 175 raise InternalError("BBOX not defined")

176

177 -class Trim(Subsetting):

178 - def __init__(self, dimension, crs, trim_low, trim_high):

179 super(Trim, self).__init__(dimension, crs) 180 181 self.trim_low = self.normalize(dimension, trim_low) 182 self.trim_high = self.normalize(dimension, trim_high)

183

184 - def validate(self, grid=None):

185 if self.trim_low is not None and self.trim_high is not None and self.trim_high < self.trim_low: 186 raise InvalidSubsettingException("Lower bound of trim greater than upper bound") 187 188 #if grid is not None: 189 # if self.dimension not in grid.axis_labels: 190 # raise InvalidAxisLabelException("Invalid axis label '%s'" % self.dimension) 191 192 if self.dimension not in ("phenomenonTime", "time", "t", "Long", "long", "Lat", "lat"): 193 raise InvalidAxisLabelException("Invalid axis label '%s'. Use 'phenomenonTime', 'Long' and 'Lat'." % self.dimension) 194 195 return True

196

197 - def getPolygon(self, footprint):

198 srid, env_minx, env_miny, env_maxx, env_maxy = self._getDataFromFootprint(footprint) 199 200 EPSILON = 1e-10 201 202 if self.dimension in ("long", "Long"): 203 miny = env_miny 204 maxy = env_maxy 205 206 if self.trim_low is None: 207 minx = env_minx 208 else: 209 minx = max(env_minx, self.trim_low) 210 211 if self.trim_high is None: 212 maxx = env_maxx 213 else: 214 maxx = min(env_maxx, self.trim_high) 215 216 elif self.dimension in ("lat", "Lat"): 217 minx = env_minx 218 maxx = env_maxx 219 220 if self.trim_low is None: 221 miny = env_miny 222 else: 223 miny = max(env_miny, self.trim_low) 224 225 if self.trim_high is None: 226 maxy = env_maxy 227 else: 228 maxy = min(env_maxy, self.trim_high) 229 230 if maxx <= minx or maxy <= miny: 231 return GEOSGeometry("POLYGON EMPTY", srid=srid) 232 else: 233 # in order to be prepared for rounding and string conversion 234 # errors, make the extent a little bit larger 235 minx = minx * (1 - EPSILON) 236 miny = miny * (1 - EPSILON) 237 maxx = maxx * (1 + EPSILON) 238 maxy = maxy * (1 + EPSILON) 239 240 poly = Polygon.from_bbox((minx, miny, maxx, maxy)) 241 poly.srid = srid 242 243 return poly

244

245 - def overlaps(self, footprint):

246 poly = self.getPolygon(footprint) 247 248 return footprint.intersects(poly)

249

250 - def contains(self, footprint):

251 poly = self.getPolygon(footprint) 252 253 return poly.contains(footprint)

254 255 #def getGridFromFile(filename, collection_srid=None): 256 # #logger.debug("CoverageInterface._getGridFromFile: SRID: %s" % str(srid)) 257 258 # return RectifiedGrid( 259 # dim=2, 260 # low=(0, 0), 261 # high=(ds.RasterXSize - 1, ds.RasterYSize - 1), 262 # axis_labels=axis_labels, 263 # srid=srid, 264 # origin=(gt[0], gt[3]), 265 # offsets=((gt[1], gt[4]), (gt[2], gt[5])) 266 # ) 267

Source Code for Module eoxserver.resources.coverages.domainset