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