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.
13 from sqlalchemy import Column
14 from sqlalchemy import Integer
15 from sqlalchemy import MetaData
16 from sqlalchemy import Table
17 from sqlalchemy import func
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 # Tests imports
25 from tests import select
26
27 metadata = MetaData()
28 Base = declarative_base(metadata=metadata)
29
30 table = Table(
31 "raster_table",
32 metadata,
33 Column("id", Integer, primary_key=True),
34 Column("geom", Geometry("POLYGON", 4326)),
35 Column("rast", Raster()),
36 )
37
38
39 class RasterTable(Base): # type: ignore
40 __tablename__ = "raster_table_orm"
41 id = Column(Integer, primary_key=True)
42 geom = Column(Geometry("POLYGON", 4326))
43 rast = Column(Raster())
44
45 def __init__(self, rast):
46 self.rast = rast
47
48
49 def test_transform_core():
50 # Define the transform query for both the geometry and the raster in a naive way
51 wrong_query = select(
52 [func.ST_Transform(table.c.geom, 2154), func.ST_Transform(table.c.rast, 2154)]
53 )
54
55 # Check the query
56 assert str(wrong_query) == (
57 "SELECT "
58 "ST_AsEWKB("
59 'ST_Transform(raster_table.geom, :ST_Transform_2)) AS "ST_Transform_1", '
60 "ST_AsEWKB(" # <= Note that the raster is processed as a Geometry here
61 'ST_Transform(raster_table.rast, :ST_Transform_4)) AS "ST_Transform_3" \n'
62 "FROM raster_table"
63 )
64
65 # Define the transform query for both the geometry and the raster in the correct way
66 correct_query = select(
67 [
68 func.ST_Transform(table.c.geom, 2154),
69 func.ST_Transform(table.c.rast, 2154, type_=Raster),
70 ]
71 )
72
73 # Check the query
74 assert str(correct_query) == (
75 "SELECT "
76 "ST_AsEWKB("
77 'ST_Transform(raster_table.geom, :ST_Transform_2)) AS "ST_Transform_1", '
78 "raster(" # <= This time the raster is correctly processed as a Raster
79 'ST_Transform(raster_table.rast, :ST_Transform_4)) AS "ST_Transform_3" \n'
80 "FROM raster_table"
81 )
82
83
84 def test_transform_ORM():
85 # Define the transform query for both the geometry and the raster in a naive way
86 wrong_query = Query([RasterTable.geom.ST_Transform(2154), RasterTable.rast.ST_Transform(2154)])
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 [
101 RasterTable.geom.ST_Transform(2154),
102 RasterTable.rast.ST_Transform(2154, type_=Raster),
103 ]
104 )
105
106 # Check the query
107 assert str(correct_query) == (
108 "SELECT "
109 "ST_AsEWKB("
110 'ST_Transform(raster_table_orm.geom, :ST_Transform_2)) AS "ST_Transform_1", '
111 "raster(" # <= This time the raster is correctly processed as a Raster
112 'ST_Transform(raster_table_orm.rast, :ST_Transform_4)) AS "ST_Transform_3" \n'
113 "FROM raster_table_orm"
114 )
