'Re: [postgis-users] Rasters in PostGIS'

[prev in list] [next in list] [prev in thread] [next in thread] 

List: postgis-users
Subject: Re: [postgis-users] Rasters in PostGIS
From: Roger_André <randre () gmail ! com>
Date: 2009-04-17 19:03:00
Message-ID: 9c2015090904171203t2a8b1736gab218becd92fb315 () mail ! gmail ! com
[Download RAW message or body]

[Attachment #2 (multipart/alternative)]

I had a project where I needed to read the pixel values from a large mosaic
for several hundred points.  The fastest method I was able to come up with
for the query was to write a Python script that did the following:

- open the raster with gdal
- use the point coordinates as input parameters to a function that
calculates the correct raster row and column based on the dimensions of the
raster, and its pixel dimensions

The IO for this tool was crude, as I used a delimited text file to contain
all of the points of interest, but it worked.

Before that, I had also "vectorized" a fairly coarse resolution raster and
stored the pixel geometries as polygons in PostGIS.  It worked very well for
doing relational queries against other data sets, but was murderously slow
for the renderer to style and display.

In both of these examples, I have not come up with a good way to extract a
specific sub-region of the data based on pixel values and return a raster.
Pulling out bounding boxes via the "gdal_translate -projwin" function is
trivial though, so maybe if all you want to do is retrieve an area that
surrounds a point, or contains a group of points, you could do that?

A combination of these methods (geometry in PostGIS, raster on disk) might
be a good solution.  I've attached the Python scripts I used for the
vectorizing, and the raster lookup, in case they might be useful starting
points.

Best of luck,

Roger
--

On Fri, Apr 17, 2009 at 8:44 AM, yoda2nd <yoda2nd@niacc.org> wrote:

> Actually I saw that WKT Raster was avalible, but I was unsure about it with
> it only in the beta stages and there does not appear to be much in the way
> of documentation yet.  I did check out the read me file in the SVN
> repository, and it seems to conflict with the documentation.
> http://trac.osgeo.org/postgis/wiki/WKTRaster/Documentation01
> https://svn.osgeo.org/postgis/spike/wktraster/README
>
> The read me file says that it will not work with the current version of
> PostGIS, but the online documentation says it will.  Not sure which one is
> correct.
>
>
>
> ------------------------------------------------------------------------
>
> I've done this before using a similar approach you have and using GDAL to
> read the elevations from the GeoTiff and dump them into an elevations table
> with point geoms.  If I had to do it again, now that we have WKT Raster in
> the works, I'd probably try that if you don't mind compiling yourself.
>
>
> http://trac.osgeo.org/postgis/wiki/WKTRaster
>
> Mateusz has a nice test drive example of WKTRaster
> http://mateusz.loskot.net/2009/03/30/crunching-wkt-rasters/
>
>
> Hope that helps,
> Regina
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users@postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

[Attachment #5 (text/html)]

I had a project where I needed to read the pixel values from a large mosaic for \
several hundred points.  The fastest method I was able to come up with for the query \
                was to write a Python script that did the following:<br><br>
- open the raster with gdal<br>
- use the point coordinates as input parameters to a function that calculates the \
correct raster row and column based on the dimensions of the raster, and its pixel \
dimensions<br><br>The IO for this tool was crude, as I used a delimited text file to \
contain all of the points of interest, but it worked.<br>

Before that, I had also &quot;vectorized&quot; a fairly coarse resolution raster \
and stored the pixel geometries as polygons in PostGIS. It worked very well for \
doing relational queries against other data sets, but was murderously slow for the \
renderer to style and display.

In both of these examples, I have not come up with a good way to extract a \
specific sub-region of the data based on pixel values and return a raster. Pulling \
out bounding boxes via the &quot;gdal_translate -projwin&quot; function is trivial \
though, so maybe if all you want to do is retrieve an area that surrounds a point, or \
contains a group of points, you could do that?

A combination of these methods (geometry in PostGIS, raster on disk) might be a \
good solution. I&#39;ve attached the Python scripts I used for the vectorizing, and \
the raster lookup, in case they might be useful starting points. Best of \
luck, Roger -- <div class="gmail_quote">On Fri, Apr 17, 2009 at 8:44 \
AM, yoda2nd &lt;<a href="mailto:yoda2nd@niacc.org" \
target="_blank">yoda2nd@niacc.org</a>&gt; wrote:

<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); \
margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Actually I saw that WKT Raster was \
avalible, but I was unsure about it with it only in the beta stages and there does \
not appear to be much in the way of documentation yet. I did check out the read me \
file in the SVN repository, and it seems to conflict with the documentation.

<a href="http://trac.osgeo.org/postgis/wiki/WKTRaster/Documentation01" \
target="_blank">http://trac.osgeo.org/postgis/wiki/WKTRaster/Documentation01</a><br> \
<a href="https://svn.osgeo.org/postgis/spike/wktraster/README" \
target="_blank">https://svn.osgeo.org/postgis/spike/wktraster/README</a><br> <br>
The read me file says that it will not work with the current version of PostGIS, but \
the online documentation says it will.  Not sure which one is correct.<br> <br>
<br>
<br>
------------------------------------------------------------------------<br>
<br>
I&#39;ve done this before using a similar approach you have and using GDAL to<br>
read the elevations from the GeoTiff and dump them into an elevations table<br>
with point geoms.  If I had to do it again, now that we have WKT Raster in<br>
the works, I&#39;d probably try that if you don&#39;t mind compiling yourself.<br>
<br>
<br>
<a href="http://trac.osgeo.org/postgis/wiki/WKTRaster" \
target="_blank">http://trac.osgeo.org/postgis/wiki/WKTRaster</a><br> <br>
Mateusz has a nice test drive example of WKTRaster <br>
<a href="http://mateusz.loskot.net/2009/03/30/crunching-wkt-rasters/" \
target="_blank">http://mateusz.loskot.net/2009/03/30/crunching-wkt-rasters/</a><br> \
<br> <br>
Hope that helps,<br>
Regina<br>
<br>
<br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@postgis.refractions.net" \
target="_blank">postgis-users@postgis.refractions.net</a><br> <a \
href="http://postgis.refractions.net/mailman/listinfo/postgis-users" \
target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br> \
</blockquote></div><br>

--0016364c716b5e98360467c4d428--

["img2shp_v2.py" (text/x-python)]

#! /usr/bin/python

# READS A RASTER AND CREATES A SHAPEFILE WITH POLYGONS 
# FOR EVERY PIXEL THAT MATCHES THE VALUE SEARCHED FOR
# IN LINE 58

import ogr
import gdal
import struct
from gdalconst import *

# POINT GDAL AT OUR INPUT DATA
filename = "sample_area.tif"
dataset = gdal.Open( filename, GA_ReadOnly )
geotransform = dataset.GetGeoTransform()

# Create x_list and y_list arrays
start_x = geotransform[0]
start_y = geotransform[3]

x_list = []
y_list = []

x_i = 0
x_val = start_x
while x_i <= dataset.RasterXSize:
  x_list.append(x_val)
  x_val += geotransform[1]
  x_i += 1

y_i = 0
y_val = start_y
while y_i <= dataset.RasterYSize:
  y_list.append(y_val)
  y_val += geotransform[5]
  y_i += 1

# CREATE A SHAPEFILE TO WRITE THE POLYS TO
driver = ogr.GetDriverByName('ESRI Shapefile')
datasource = driver.CreateDataSource('sample_area.shp')
layer = datasource.CreateLayer('polys', geom_type=ogr.wkbPolygon)

# CREATE 'PIXEL_VALUE' ATTRIBUTE FIELD
id_field = ogr.FieldDefn('PIXEL_VALUE', ogr.OFTString)	
layer.CreateField(id_field)						

# READ RASTER DATA
band = dataset.GetRasterBand(1)
row = 0
polys = 0
while row < dataset.RasterYSize:
  scanline = band.ReadRaster( 0, row, band.XSize, 1, band.XSize, 1, GDT_Float32 )
  tuple_of_floats = struct.unpack('f' * band.XSize, scanline)
  total_cols = len(tuple_of_floats)
  col = 0
  for col in range(0,total_cols):
    pixel = tuple_of_floats[col]
    if pixel != 0.0: # <-------------- Set your pixel value here.
      print 'Creating poly for Col: %d, Row: %d' % (col, row)
      xi = col
      yi = row

      # BUILD POLYGON GEOMETRY USING VALUES FROM x_list AND y_list ARRAYS
      x1 = x_list[xi]
      x2 = x_list[xi]
      x3 = x_list[xi + 1]
      x4 = x_list[xi + 1]

      y1 = y_list[yi]
      y2 = y_list[yi + 1]
      y3 = y_list[yi + 1]
      y4 = y_list[yi]
      wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))' % (x1,y1, x2,y2, x3,y3, x4,y4, x1,y1)

      # THIS IS WHERE WE CREATE A POLYGON FOR THAT PIXEL
      geom = ogr.CreateGeometryFromWkt(wkt)
      feat = ogr.Feature(layer.GetLayerDefn())
      feat.SetGeometry(geom)

      # SET THE 'PIXEL_VALUE' ATTRIBUTE VALUE
      feat.SetField('PIXEL_VALUE', pixel) 

      # CREATE THE NEW FEATURE
      layer.CreateFeature(feat)
      polys += 1
  row += 1

# CLEAN UP AFTER SHAPEFILE CREATION
datasource.Destroy()

print "Created %d polys" % (polys)

["raster_intersect.py" (text/x-python)]

#! /usr/bin/python

"""Takes a point and calculates what the pixel value is of an intersecting raster."""

# USAGE: r.intersect.py <raster> <points_csv> <new_field_name>
# ASSUMES DATA IS IN THIS FORMAT
# pop|name|state|fips|lon lat

import sys
import numpy
from osgeo import gdal

def testExtents(point, ul, lr):
  # TESTS TO SEE IF POINT IS WITHIN THE EXTENTS OF THE RASTER
  point_x, point_y = point
  x_min, y_max = ul
  x_max, y_min = lr
  if not x_min <= point_x <= x_max or not y_min <= point_y <= y_max:
    return 0
  else:
    return 1

def calcIndex(point, pix_x, ul_x, pix_y, ul_y):
  # CALCULATE ARRAY POSITION OF POINT IN RASTER
  point_x, point_y = point
  x_diff = point_x - ul_x
  y_diff = ul_y - point_y
  x_pixels = abs(x_diff/pix_x)
  y_pixels = abs(y_diff/pix_y)
  x_off = int(x_pixels)
  y_off = int(y_pixels)
  offset = [x_off, y_off]
  return offset

def main():
  src_file = sys.argv[1]
  cities_file = sys.argv[2]
  new_field = sys.argv[3]

  # GATHER BASIC INFO ABOUT RASTER
  src_ds = gdal.Open( src_file )
  numbands = src_ds.RasterCount
  if numbands > 1:
    print "ERROR: Raster must be single-band"
    sys.exit(1)
  src_band = src_ds.GetRasterBand(1)
  src_geotransform = src_ds.GetGeoTransform()
  width = src_ds.RasterXSize
  pix_x = float(src_geotransform[1])
  ul_x = float(src_geotransform[0])
  height = src_ds.RasterYSize
  pix_y = float(src_geotransform[5])
  ul_y = float(src_geotransform[3])
  lr_x = (pix_x * width) + ul_x
  lr_y = (pix_y * height) + ul_y
  ul = [ul_x, ul_y]
  lr = [lr_x, lr_y]

  # READ POINT DATA, ASSUMES DATA IS IN THIS FORMAT
  # pop|name|state|fips|lon lat
  city_data = open(cities_file).readlines()
  print "%s|%s" % (city_data[0].strip(), new_field)
  for city in city_data[1:]:
    data_elem = city.split("|")
    lon_lat = data_elem[4].split()
    lon, lat = lon_lat
    point = (float(lon), float(lat))
    if testExtents(point, ul, lr):
      offset = calcIndex(point, pix_x, ul_x, pix_y, ul_y)
      # GET RASTER VALUE AT SPECIFIC ROW/COLUMN OFFSET IN RASTER
      iX = offset[0]
      iY = offset[1]
      src_data = src_band.ReadAsArray(0,iY,src_ds.RasterXSize,1)
      col_values = src_data[0]
      raster_value = col_values[iX]
      print "%s|%f" % (city.strip(), raster_value)
    else:
      print "%s|no_data" % (city.strip())

if __name__ == "__main__":

  main()

_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[prev in list] [next in list] [prev in thread] [next in thread] 