Changeset 3920 for spike/wktraster/scripts/gdal2wktraster.py

Timestamp:

03/21/09 09:55:58 (3 years ago)

Author:

mloskot

Message:

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.

Files:

• spike/wktraster/scripts/gdal2wktraster.py (modified) (3 diffs)

spike/wktraster/scripts/gdal2wktraster.py

@@ -274 +274 @@
     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):
@@ -437 +460 @@
     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)
@@ -573 +596 @@
         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],