Source Code for Module madrona.raster_stats.models

  1  from django.db import models 
  2  from django.conf import settings 
  3  from django.core import serializers 
  4  from django.db.utils import DatabaseError 
  5  from madrona.common.utils import get_logger 
  6  import tempfile 
  7  import time 
  8  import os 
  9   
 10  logger = get_logger() 
 11   
 12  try: 
 13      RASTDIR = settings.RASTER_DIR 
 14  except: 
 15      RASTDIR = os.path.join(os.path.dirname(__file__), 'test_data') 
 16   
 17  RASTER_TYPES = (  
 18                  ("continuous", "continuous"), 
 19                  ("categorical", "catgorical"), 
 20                 ) 
 21   
 22  STARSPAN_BIN = settings.STARSPAN_BIN 
 23   
 24 -class RasterDataset(models.Model): 
 25      name = models.CharField(max_length=30, unique=True) 
 26      full_name = models.CharField(max_length=255, default="") 
 27      filepath = models.FilePathField(max_length=255, path=RASTDIR, recursive=True) 
 28      type = models.CharField(max_length=30, choices=RASTER_TYPES) 
 29   
 30 -    def __unicode__(self): 
 31          return unicode(self.name + " raster at " + self.filepath) 
 32   
 33 -    def save(self, *args, **kwargs): 
 34          super(RasterDataset, self).save(*args, **kwargs) 
 35          # invalidate the zonal stats cache 
 36          caches = ZonalStatsCache.objects.filter(raster=self) 
 37          caches.delete() 
 38   
 39      @property 
 40 -    def is_valid(self): 
 41          # TODO is there a CHEAP way to check for GDAL open? 
 42          if not os.path.exists(self.filepath): 
 43              return False 
 44          return True 
 45   
 46 -class ZonalCategory(models.Model): 
 47      category = models.IntegerField() 
 48      count = models.IntegerField() 
 49   
 50 -    def __unicode__(self): 
 51          return unicode("Category %s (count = %s)" % (self.category, self.count)) 
 52   
 53 -class ZonalStatsCache(models.Model): 
 54      geom_hash = models.CharField(max_length=255) 
 55      raster = models.ForeignKey('RasterDataset') 
 56      sum = models.FloatField(null=True, blank=True) 
 57      avg = models.FloatField(null=True, blank=True) 
 58      min = models.FloatField(null=True, blank=True) 
 59      max = models.FloatField(null=True, blank=True) 
 60      mode = models.FloatField(null=True, blank=True) 
 61      median = models.FloatField(null=True, blank=True) 
 62      stdev = models.FloatField(null=True, blank=True) 
 63      nulls = models.FloatField(null=True, blank=True) 
 64      pixels = models.FloatField(null=True, blank=True) 
 65      categories = models.ManyToManyField(ZonalCategory) 
 66      date_modified = models.DateTimeField(auto_now=True) 
 67   
 68      @property 
 69 -    def json(self): 
 70          return serializers.serialize("json", self) 
 71   
 72 -    def __unicode__(self): 
 73          return unicode("Zonal Stats for %s - avg:%s , pixels:%s, nulls:%s" % (self.raster.name, self.avg, self.pixels, self.nulls)) 
 74   
 75 -    class Meta: 
 76          unique_together = ('geom_hash', 'raster') 
 77   
 78 -def geom_to_file(geom, filepath): 
 79      json = """{ 
 80  "type": "FeatureCollection", 
 81  "features": [ 
 82  { "type": "Feature", "properties": { "id": 1 }, "geometry": %s } 
 83  ] 
 84  }""" % geom.json 
 85      fh = open(filepath,'w') 
 86      fh.write(json) 
 87      fh.close() 
 88      assert os.path.exists(filepath) 
 89   
 90 -def _run_starspan_zonal(geom, rasterds, pixprop=0.5): 
 91      """ 
 92      Consider this a 'private' method .. dont call directly, use zonal_stats() instead 
 93      Runs starspan and returns a ZonalStatsCache object 
 94      """ 
 95      # Create tempdir and cd in  
 96      tmpdir_base = tempfile.gettempdir() 
 97      geom_hash = geom.wkt.__hash__() 
 98      timestamp = str(time.time()) 
 99      tmpdir = os.path.join(tmpdir_base, 'madrona.raster_stats', str(geom_hash), timestamp, rasterds.name) 
100      os.makedirs(tmpdir) 
101      old_dir = os.getcwd() 
102      os.chdir(tmpdir) 
103   
104      # Output geom to temp dataset 
105      out_json = os.path.join(tmpdir, 'geom_%s.json' % timestamp) 
106      geom_to_file(geom, out_json) 
107   
108      # Run starspan 
109      out_csv = os.path.join(tmpdir, 'output_%s_stats.csv' % timestamp) 
110      out_categories = os.path.join(tmpdir, 'output_%s_categories.csv' % timestamp) 
111      if os.path.exists(out_csv):  
112          os.remove(out_csv) 
113      if os.path.exists(out_categories):  
114          os.remove(out_categories) 
115      cmd = '%s --vector %s --where "id=1" --out-prefix %s/output_%s --out-type table --summary-suffix _stats.csv --raster %s --stats avg mode median min max sum stdev nulls --pixprop %s ' % (STARSPAN_BIN,out_json,tmpdir, timestamp, rasterds.filepath, pixprop) 
116      if rasterds.type == 'categorical': 
117          cmd += "--class-summary-suffix _categories.csv" 
118   
119      starspan_out = os.popen(cmd).read() 
120   
121      if not os.path.exists(out_csv): 
122          raise Exception("Starspan failed to produce output file: %s" % starspan_out) 
123      if rasterds.type == 'categorical' and not os.path.exists(out_categories): 
124          raise Exception("Starspan failed to produce output file: %s" % starspan_out) 
125   
126      res = open(out_csv,'r').readlines() 
127   
128      # Create zonal model 
129      zonal, created = ZonalStatsCache.objects.get_or_create(raster=rasterds, geom_hash=geom_hash) 
130   
131      # Make sure we have valid results output by starspan 
132      if len(res) == 2 and "Intersecting features: 0" not in starspan_out: 
133          headers = [x.strip() for x in res[0].split(',')] 
134          vals = [x.strip() for x in res[1].split(',')] 
135          assert len(headers) == len(vals) 
136   
137          # loop and populate model 
138          for i in range(len(headers)): 
139              if "_Band1" in headers[i]: 
140                  stat_type = headers[i].replace("_Band1",'') 
141                  zonal.__dict__[stat_type] = float(vals[i]) 
142              elif headers[i] == 'numPixels': 
143                  zonal.pixels = float(vals[i]) 
144   
145      # Handle categories 
146      if rasterds.type == 'categorical' and zonal.pixels: 
147          zonal.save() 
148          res = open(out_categories).readlines() 
149          headers = [x.strip() for x in res[0].split(',')] 
150          assert headers == ['FID', 'class', 'count'] 
151          for row in res[1:]: 
152              vals = [int(x.strip()) for x in row.split(',')] 
153              cat = vals[1] 
154              count = vals[2] 
155              zc = ZonalCategory.objects.create(category=cat, count=count) 
156              zc.save() 
157              zonal.categories.add(zc) 
158   
159      try: 
160          if zonal.pixels: 
161              zonal.save() 
162      except: 
163          # Most likely another zonal stats cache for this geom/raster 
164          # was saved to the cache before this one completed. 
165          pass 
166   
167      if settings.STARSPAN_REMOVE_TMP: 
168          os.chdir(old_dir) 
169          import shutil 
170          shutil.rmtree(tmpdir) 
171   
172      return zonal 
173   
174 -def clear_cache(): 
175      objs = ZonalStatsCache.objects.all() 
176      objs.delete() 
177   
178 -def zonal_stats(geom, rasterds, read_cache=True, cache_only=False, pixprop=0.5): 
179      """ 
180      Given a GEOSGeometry and a RasterDataset, 
181      compute the zonal stats and return json like 
182       { 'raster': 'elevation', 'stats': {'sum': 10234.2, 'mean': 12.4}} 
183      result will be stored in cache 
184       and cache value is returned if read_cache 
185      """ 
186      if not geom.valid: 
187          return None 
188   
189      geom_hash = geom.wkt.__hash__() 
190   
191      cached = None 
192      if read_cache: 
193          try: 
194              cached = ZonalStatsCache.objects.get(geom_hash=geom_hash, raster=rasterds) 
195          except ZonalStatsCache.DoesNotExist: 
196              cached = None 
197          except DatabaseError: 
198              cached = None 
199   
200      if cached: 
201          result = cached 
202          result.from_cache = True 
203      else: 
204          if cache_only: 
205              # Return an empty result 
206              result = ZonalStatsCache(geom_hash=geom_hash, raster=rasterds) 
207          else: 
208              result = _run_starspan_zonal(geom, rasterds, pixprop=pixprop) 
209          result.from_cache = False 
210   
211      return result 
212