[postgis-commits] svn - r3920 - spike/wktraster/scripts

[postgis-commits at postgis.refractions.net]

Author: mloskot
Date: 2009-03-21 09:55:58 -0700 (Sat, 21 Mar 2009)
New Revision: 3920

Modified:
   spike/wktraster/scripts/gdal2wktraster.py
Log:
gdal2wktraster: Calculate real raster extent for regular blocking and insert it to RASTER_COLUMNS. Added calculate_geoxy() and calculate_bounding_box() functions. This is much more usable/stable version of raster blocking support. All loading modes tested using GDAL test dataset: utm.tif, 512x512, 1 band, 1BUI. TODO: Test with multi-band datasets.

Modified: spike/wktraster/scripts/gdal2wktraster.py
===================================================================
--- spike/wktraster/scripts/gdal2wktraster.py	2009-03-21 15:59:56 UTC (rev 3919)
+++ spike/wktraster/scripts/gdal2wktraster.py	2009-03-21 16:55:58 UTC (rev 3920)
@@ -273,7 +273,30 @@
     ny = raster_size[1] / block_size[1]
     return (int(nx), int(ny))
 
+def calculate_geoxy(ds, xy):
+    """Calculate georeferenced coordinate from given x and y"""
+    assert ds is not None
+    assert xy is not None
+    assert len(xy) == 2
 
+    gt = ds.GetGeoTransform()
+    xgeo = gt[0] + gt[1] * xy[0] + gt[2] * xy[1];
+    ygeo = gt[3] + gt[4] * xy[0] + gt[5] * xy[1];
+
+    return (xgeo, ygeo)
+
+def calculate_bounding_box(ds):
+    """Calculate georeferenced coordinates of spatial extent of raster dataset"""
+    assert ds is not None
+
+    # UL, LL, UR, LR
+    dim = ( (0,0),(0,ds.RasterYSize),(ds.RasterXSize,0),(ds.RasterXSize,ds.RasterYSize) )
+
+    ext = (calculate_geoxy(ds, dim[0]), calculate_geoxy(ds, dim[1]),
+           calculate_geoxy(ds, dim[2]), calculate_geoxy(ds, dim[3]))
+
+    return ext
+
 def check_hex(hex, bytes_size = None):
     assert hex is not None, "Error: Missing hex string"
     size = len(hex)
@@ -436,7 +459,7 @@
     pixel_size = ( ds.GetGeoTransform()[1], ds.GetGeoTransform()[1] )
     pixel_types = collect_pixel_types(ds, band_from, band_to)
     nodata_values = collect_nodata_values(ds, band_from, band_to)
-    extent = ((0,0), (1,0), (1,1), (0,1)) # TODO: Generate polygon geometry containing raster/all tiles
+    extent = calculate_bounding_box(ds)
     sql = make_sql_addrastercolumn(options, pixel_types, nodata_values, pixel_size, block_size, extent)
     options.output.write(sql)
 
@@ -572,7 +595,7 @@
         assert len(extent[0]) == len(extent[3]) == 2, "Invalid extent, pair of X and Y pair expected"
         rb = 'true'
         bs = ( blocksize[0], blocksize[1] )
-        extgeom = "ST_GeometryFromText('POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))', %d)" % \
+        extgeom = "ST_Box2D(ST_GeometryFromText('POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))', %d))" % \
                   (extent[0][0], extent[0][1], extent[1][0], extent[1][1],
                    extent[2][0], extent[2][1], extent[3][0], extent[3][1],
                    extent[0][0], extent[0][1], options.srid)