4

I'm trying to convert a raster to a polygon and export it from PostGIS as GeoJSON. Specifically I want each pixel to be a polygon with associated the raster value.

PostGIS seems to have a few different functions for converting rasters to polygons, I'm using the ST_PixelAsPolygons function because I want each pixel to be represented as its own polygon.

Example:

SELECT ST_PixelAsPolygon(rast) FROM tbl;

However this returns a set of records for each raster and I can't easily convert this to GeoJSON.

Does anyone know of a simple way to do this?

7

If you look at the example on the man page for ST_PixelAsPolygons you will see how you can access the geometries using table_alias.geom syntax (similar in spirit to how ST_Dump works to turn a set or records into individual rows). Following on from that example, you can pass (gv).geom to the ST_AsGeoJSON function, eg,

SELECT ST_AsGeoJSON((gv).geom) 
  FROM 
      (SELECT ST_PixelAsPolygons(
          ST_SetValue(
              ST_SetValue(
                  ST_AddBand(
                    ST_MakeEmptyRaster(2, 2, 0, 0, 0.001, -0.001, 0.001, 0.001, 4269), 
                                       '8BUI'::text, 1, 0), 
                    2, 2, 10), 
            1, 1, NULL)
     ) gv 
 ) foo;

which will return the three rows:

{"type":"Polygon","coordinates":[[[0.001,0.001],[0.002,0.002],[0.003,0.001],[0.002,0],[0.001,0.001]]]}

{"type":"Polygon","coordinates":[[[0.001,-0.001],[0.002,0],[0.003,-0.001],[0.002,-0.002],[0.001,-0.001]]]}

{"type":"Polygon","coordinates":[[[0.002,0],[0.003,0.001],[0.004,0],[0.003,-0.001],[0.002,0]]]}

If you wish to turn this into a FeatureCollection, you can use a combination of row_to_json, array_agg and array_to_json to group the rows together into an array and add in the "type":"FeatureCollection" before the features array. The general idea is that you select the 'Feature', the GeoJSON representation of the geometry and the properties in one subquery, and then use array_to_json and array_agg to group this into an array of features, with 'FeatureCollection' as type before the array. Finally, you select row_to_json to group this all together inside {}.

Note you have to do the nasty subselect on the properties -- in your case the raster pixel values -- to avoid the labels being f1, f2, etc. See this Postgres Online article for more information.

SELECT row_to_json(fc) 
  FROM 
    (SELECT 'FeatureCollection' as type, array_to_json(array_agg(feats)) as features 
      FROM 
         (SELECT 'Feature' as type, 
                 st_asgeojson((gv).geom)::json as geometry, 
                 row_to_json((SELECT props FROM (SELECT (gv).val as value) as props )) as properties  
            FROM 
                (SELECT 
                    ST_PixelAsPolygons(
                       ST_SetValue(
                         ST_SetValue(
                             ST_AddBand(
                                 ST_MakeEmptyRaster(2, 2, 0, 0, 0.001, -0.001, 0.001, 0.001, 4269), 
                            '8BUI'::text, 1, 0),
                        2, 2, 10), 
                     1, 1, NULL)
               ) gv 
          ) json
    ) feats 
) fc;

which now produces:

{"type":"FeatureCollection", 
     "features":[
           {"type":"Feature",
                "geometry":{"type":"Polygon","coordinates":[[[0.001,0.001],[0.002,0.002],[0.003,0.001],[0.002,0],[0.001,0.001]]]},
                "properties":{"value":1}
           },
           {"type":"Feature",
                "geometry":{"type":"Polygon","coordinates":[[[0.001,-0.001],[0.002,0],[0.003,-0.001],[0.002,-0.002],[0.001,-0.001]]]},
                "properties":{"value":1}
            },
            {"type":"Feature",
                "geometry":{"type":"Polygon","coordinates":[[[0.002,0],[0.003,0.001],[0.004,0],[0.003,-0.001],[0.002,0]]]},
                "properties":{"value":10}
            }
     ]
}

Note: You need to cast the ST_AsGeoJSON to JSON, with the :: operator, or else you will end up with all the quotes being escaped, as \" in the output.

  • Did I miss the point of your question? I realize that this is perhaps, not simple, but I believe it answers the question? – John Powell Jan 23 '16 at 9:53
  • OK, I did miss something. I have edited the post which now includes the pixel values as a property of the output GeoJSON. Obviously, the whole inner select making a dummy raster would be replaced with a straight call to ST_PixelAsPolygons in your case. – John Powell Jan 26 '16 at 10:00
  • I have looked at the solution you provided (thanks by the way). Here I have pasted the query I want to run. I need all the pixels from a raster that intersects a polygon in another table. The query works but I get to many polygons in the GeoJSON. I have posted an image where I've outlined the polygon queried in the nested select clause, find it here. I would be greatful for any advice, even though it's not really related to the original question. – Simtrax Oct 1 '17 at 9:30
  • OK, you should probably accept this answer and ask another question :-). However, the function ST_Reclass allows you to use banding, so that when you convert rasters to polygons you don't get so many. – John Powell Oct 1 '17 at 12:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.