checkpoint cause I'm leaving (fails tests)

geoviz-tsv

@@ -7,8 +7,13 @@ foo.geo.tsv, it emits:
    1. a GeoJSON file containing each tweet as a point (foo.geojson).
    2. a GeoTIFF with a heat map with tweet counts in each pixel (foo.tiff).'''
 
+from __future__ import division
+
 import argparse
 
+import numpy as np
+
+from geo import srs
 import testable
 import tweet
 import u
@@ -30,7 +35,91 @@ ap.add_argument('--verbose',
 ap.add_argument('input',
                 metavar='FILE',
                 help='TSV file containing geocoded tweets')
-
+
+
+### Classes ###
+
+
+class GeoTIFF(object):
+   '''Create a GeoTIFF file from a NumPy array and write it to disk. E.g.:
+   
+      >>> import tempfile
+      >>> d = np.array([[1, 0.5, 0.5], [0.5, 0.5, 0.25]])
+      >>> f = GeoTIFF(srs.NWMAX, srs.SEMAX, d)
+      >>> fn = tempfile.mktemp()  # insecure
+      >>> f.write(fn + '.tif')
+      >>> os.getsize(fn) == 1234
+      
+      Note that the above doesn't verify that the GeoTIFF is actually correct.
+      What you should see (and the above leaves the temp file laying around) is
+      a six-pixel image across the whole map with the northwest pixel darker
+      and the southeast pixel lighter. The pixels will not be square.'''
+   pass
+
+   def __init__(self, nw, se, data=None):
+      '''Parameters:
+      
+           nw ..... Northwest corner of geodata (geopoint)
+           se ..... Southwest corner (geopoint)
+           data ... Grid to write (2-D NumPy array) (optional)
+         
+         nw and se must have a valid (and the same) SRS.'''
+      if (nw.srid is None or nw.srid != se.srid):
+         raise ValueError('one or more invalid SRS: (%s, %s)'
+                          % (nw.srid, se.srid))
+      if (not srs.is_northwest(nw, se)):
+         raise ValueError('corner %s must be northwest of %s' % (nw, se))
+      self.nw = nw
+      self.se = se
+      if (data is not None):
+         self.data = data
+
+   @property
+   def data(self):
+      return self._data
+
+   @data.setter
+   def data(self, value):
+      if (len(value.shape) != 2):
+         raise ValueError('data dimensions must be 2, not %d' % (len(value.shape)))
+
+   def height_for_square(self, width):
+      '''Return the height which gives square pixels for width. (Note that
+         squareness only holds if the image is not reprojected before display.)'''
+      assert False, 'unimplemented'
+
+   def write(self, filename, type='Float32', jpeg=False):
+      'Write a GeoTIFF to filename with data type type.'
+      assert (not jpeg), 'unimplemented'
+      if (jpeg):
+         assert False, 'unimplemented'
+         if (type != 'Byte'):
+            raise ValueError('jpeg requires Byte type, not %s' % (type))
+         options = ['COMPRESS=JPEG', 'JPEG_QUALITY=95']
+      else:
+         options = ['COMPRESS=DEFLATE', 'PREDICTOR=2']
+      # Set up GDAL driver
+      driver = ogdal.GetDriverByName('GTiff')
+      out = driver.create(filename, data.shape[0], data.shape[1],
+                          1, getattr(ogdal, 'GDT_' + type), options)
+      # Affine transform from image space to projected space. I don't quite
+      # understand what is going on here; the resulting image has upper left
+      # and lower left corners reversed according to gdalinfo (and same for
+      # right). However, it displays fine in QGIS. An alternative is to offer
+      # ymax and invert the pixel size (making it negative), which gives
+      # corners that seem right but then the image is upside down.
+      # http://gdal.org/classGDALDataset.html#af9593cc241e7d140f5f3c4798a43a668
+      out.SetGeoTransform([xmin, (xmax - xmin) / width_px, 0,
+                           ymin, 0, (ymax - ymin) / height_px])
+      out.SetProjection(srs.wkt(nw.srid))
+      # Write the main data.
+      out.GetRasterBand(1).WriteArray(probs)
+      # In order to correctly display in QGIS, you need to compute the exact
+      # statistics. A bug prevents QGIS from doing this
+      # http://hub.qgis.org/issues/6496(), and also if use this call, then
+      # it's embedded in the file and no auxiliary .xml file is created.
+      out.GetRasterBand(1).GetStatistics(0,1)
+
 
 ### Main ###