1""" 2Reproject a Raster using ST_Transform 3===================================== 4 5The `ST_Transform()` function (and a few others like `ST_SnapToGrid()`) can be used on 6both `Geometry` and `Raster` types. In `GeoAlchemy2`, this function is only defined for 7`Geometry` as it can not be defined for several types at the same time. Thus using this 8function on `Raster` requires minor tweaking. 9 10This example uses both SQLAlchemy core and ORM queries. 11""" 12from sqlalchemy import Column 13from sqlalchemy import func 14from sqlalchemy import Integer 15from sqlalchemy import MetaData 16from sqlalchemy import select 17from sqlalchemy import Table 18from sqlalchemy.orm import Query 19from sqlalchemy.ext.declarative import declarative_base 20 21from geoalchemy2 import Geometry 22from geoalchemy2 import Raster 23 24 25metadata = MetaData() 26Base = declarative_base(metadata=metadata) 27 28table = Table( 29 "raster_table", 30 metadata, 31 Column("id", Integer, primary_key=True), 32 Column("geom", Geometry("POLYGON", 4326)), 33 Column("rast", Raster(srid=4326)), 34) 35 36 37class RasterTable(Base): 38 __tablename__ = 'raster_table_orm' 39 id = Column(Integer, primary_key=True) 40 geom = Column(Geometry("POLYGON", 4326)) 41 rast = Column(Raster(srid=4326)) 42 43 def __init__(self, rast): 44 self.rast = rast 45 46 47def test_transform_core(): 48 # Define the transform query for both the geometry and the raster in a naive way 49 wrong_query = select([ 50 func.ST_Transform(table.c.geom, 2154), 51 func.ST_Transform(table.c.rast, 2154) 52 ]) 53 54 # Check the query 55 assert str(wrong_query) == ( 56 "SELECT " 57 "ST_AsEWKB(" 58 "ST_Transform(raster_table.geom, :ST_Transform_2)) AS \"ST_Transform_1\", " 59 "ST_AsEWKB(" # <= Note that the raster is processed as a Geometry here 60 "ST_Transform(raster_table.rast, :ST_Transform_4)) AS \"ST_Transform_3\" \n" 61 "FROM raster_table" 62 ) 63 64 # Define the transform query for both the geometry and the raster in the correct way 65 correct_query = select([ 66 func.ST_Transform(table.c.geom, 2154), 67 func.ST_Transform(table.c.rast, 2154, type_=Raster) 68 ]) 69 70 # Check the query 71 assert str(correct_query) == ( 72 "SELECT " 73 "ST_AsEWKB(" 74 "ST_Transform(raster_table.geom, :ST_Transform_2)) AS \"ST_Transform_1\", " 75 "raster(" # <= This time the raster is correctly processed as a Raster 76 "ST_Transform(raster_table.rast, :ST_Transform_4)) AS \"ST_Transform_3\" \n" 77 "FROM raster_table" 78 ) 79 80 81def test_transform_ORM(): 82 # Define the transform query for both the geometry and the raster in a naive way 83 wrong_query = Query([ 84 RasterTable.geom.ST_Transform(2154), 85 RasterTable.rast.ST_Transform(2154) 86 ]) 87 88 # Check the query 89 assert str(wrong_query) == ( 90 "SELECT " 91 "ST_AsEWKB(" 92 "ST_Transform(raster_table_orm.geom, :ST_Transform_2)) AS \"ST_Transform_1\", " 93 "ST_AsEWKB(" # <= Note that the raster is processed as a Geometry here 94 "ST_Transform(raster_table_orm.rast, :ST_Transform_4)) AS \"ST_Transform_3\" \n" 95 "FROM raster_table_orm" 96 ) 97 98 # Define the transform query for both the geometry and the raster in the correct way 99 correct_query = Query([ 100 RasterTable.geom.ST_Transform(2154), 101 RasterTable.rast.ST_Transform(2154, type_=Raster) 102 ]) 103 104 # Check the query 105 assert str(correct_query) == ( 106 "SELECT " 107 "ST_AsEWKB(" 108 "ST_Transform(raster_table_orm.geom, :ST_Transform_2)) AS \"ST_Transform_1\", " 109 "raster(" # <= This time the raster is correctly processed as a Raster 110 "ST_Transform(raster_table_orm.rast, :ST_Transform_4)) AS \"ST_Transform_3\" \n" 111 "FROM raster_table_orm" 112 ) 113