How to properly rotate a raster

Sat Mar 23 02:04:48 PDT 2024

Found it:

gdal_translate -of GTiff PG:"host=xxxxxxx port=5432 dbname=xxxxxxx user=xxxxxxx schema=xxxxxxx table=rast_not_rotated mode=2" /home/xxxxxxx/Downloads/rast_not_rotated.tif
Input file size is 1000, 1000
0...10...20...30...40...50...60...70...80...90...100 - done.

gdal_translate -of GTiff PG:"host=xxxxxxx port=5432 dbname=xxxxxxx user=xxxxxxx schema=xxxxxxx table=rast_rotated_from_raster mode=2" /home/xxxxxxx/Downloads/rast_rotated_from_raster.tif
ERROR 1: GDAL PostGIS Raster driver can not work with rotated rasters yet.
Input file size is -2147483648, -2147483648
ERROR 1: Error: -srcwin 0 0 -2.14748e+09 -2.14748e+09 has negative width and/or height.

gdal_translate -of GTiff PG:"host=xxxxxxx port=5432 dbname=xxxxxxx user=xxxxxxx schema=xxxxxxx table=rast_rotated_from_vector mode=2" /home/xxxxxxx/Downloads/rast_rotated_from_vector.tif
Input file size is 1000, 1000
0...10...20...30...40...50...60...70...80...90...100 - done. All good here, raster nicely rotated.

So, I better rotate the vector layers first and then convert them into rasters. At least until rotated rasters are implemented in GDAL driver.

And the reason why DBeaver was showing it nicely rotated was because we were looking at a Polygon created from the rotated raster and PostGIS Raster does know how to deal with rotated rasters.

Good to have this cleared out.

Cheers!

Eloi

On Friday, March 22nd, 2024 at 20:06, Regina Obe <lr at pcorp.us> wrote:

> The only reason I can think why QGIS and DBeaver would show something different is how they are projecting.
>
> I’m guessing they are both projecting to Web Mercator or some such thing to overlay your geometry on the map, and their assumptions about SRID 28992 are different.
>
> I’d try transforming both your geometries to 4326
>
> ST_Transform( ST_Polygon(…), 4326)
>
> I’m guessing whatever spatial projecting they do from 4326 should be consistent.
>
> Alternatively ditch the map overlay and then maybe they won’t try to reproject.
>
> From: Eloi Ribeiro <mail at eloiribeiro.eu>
> Sent: Thursday, March 21, 2024 11:47 AM
> To: PostGIS Users Discussion <postgis-users at lists.osgeo.org>
> Subject: RE: How to properly rotate a raster
>
> Hi Regina and all,
>
> I' seeing different things depending on what application I use. On QGIS the rotated raster appears smaller in width and height. On DBeaver, both rotate, although in different directions.
>
> -- Create raster from polygon (black line)
>
> DROP TABLE IF EXISTS public.rast_not_rotated;
>
> CREATE TABLE public.rast_not_rotated AS
>
> SELECT ST_AsRaster(t.geom,1000,1000,'8BUI',1,0) AS rast
>
> FROM (SELECT ST_GeomFromText('POLYGON((166989.36462303242 505764.77048559673,
>
> 166996.25788353998 505576.476318614,
>
> 167265.54134611748 505585.94966528506,
>
> 167258.67088940175 505774.2438322678,
>
> 166989.36462303242 505764.77048559673))',28992) geom) t;
>
> ALTER TABLE public.rast_not_rotated ADD COLUMN rid SERIAL PRIMARY KEY;
>
> CREATE INDEX rast_not_rotated_rast_gist ON public.rast_not_rotated USING gist (st_convexhull(rast));
>
> SELECT AddRasterConstraints('public'::name, 'rast_not_rotated'::name,'rast'::name);
>
> -- rotate previous raster (blue)
>
> DROP TABLE IF EXISTS public.rast_rotated_from_raster;
>
> CREATE TABLE public.rast_rotated_from_raster AS SELECT ST_SetRotation(rast, -0.5) rast FROM public.rast_not_rotated;
>
> ALTER TABLE public.rast_rotated_from_raster ADD COLUMN rid SERIAL PRIMARY KEY;
>
> CREATE INDEX rast_rotated_from_raster_rast_gist ON public.rast_rotated_from_raster USING gist (st_convexhull(rast));
>
> SELECT AddRasterConstraints('public'::name, 'rast_rotated_from_raster'::name,'rast'::name);
>
> -- rotate polygon and then save as raster (red)
>
> DROP TABLE IF EXISTS public.rast_rotated_from_vector;
>
> CREATE TABLE public.rast_rotated_from_vector AS
>
> SELECT ST_AsRaster(t.geom,1000,1000,'8BUI',1,0) AS rast
>
> FROM (SELECT ST_Rotate(
>
> ST_GeomFromText('POLYGON((166989.36462303242 505764.77048559673,
>
> 166996.25788353998 505576.476318614,
>
> 167265.54134611748 505585.94966528506,
>
> 167258.67088940175 505774.2438322678,
>
> 166989.36462303242 505764.77048559673))',28992)
>
> ,-0.5, 166989.36462303242, 505764.77048559673) geom) t;
>
> ALTER TABLE public.rast_rotated_from_vector ADD COLUMN rid SERIAL PRIMARY KEY;
>
> CREATE INDEX rast_rotated_from_vector_rast_gist ON public.rast_rotated_from_vector USING gist (st_convexhull(rast));
>
> SELECT AddRasterConstraints('public'::name, 'rast_rotated_from_vector'::name,'rast'::name);
>
> On QGIS.
>
> black line - raster from polygon
>
> red - the rotated polygon and then converted to raster
>
> blue - rotated raster, different from how it is displayed in DBeaver.
>
> On DBeaver this is how it shows the rotated raster.
>
> On DBeaver this is how it shows the rotated polygon and then converted to raster.
>
> On Tuesday, March 19th, 2024 at 14:58, Regina Obe <lr at pcorp.us> wrote:
>
>> An easier way to compare what is happening is to look at the ST_Polygon output of the two. https://postgis.net/docs/en/RT_ST_Polygon.html
>>
>> I’m thinking the confusion (which I had thought myself) is that the width and height is measured against the X y axis, but now I think about it, I think it’s measured along the skew so width and height would remain the same in both cases. In the case of your polygon you rotate and then save as raster, the new raster is axis aligned and so there you do see the width and height change.
>>
>> So when I run this query:
>>
>> SELECT ST_Polygon(r1.rast) AS orig_geom_rast, ST_Polygon(r2.rast) AS rotated_geom_rast, ST_Width(r1.rast) AS orig_width, ST_Width(r2.rast) AS rotated_width
>>
>> FROM (SELECT rast FROM rast_not_rotated) AS r1,
>>
>> (SELECT ST_SetRotation(rast, -0.5 ) AS rast FROM rast_not_rotated) AS r2;
>>
>> The dashed polygon is the original polygon and the filled in one is after you rotate the raster. In both cases the ST_Width and ST_Height do not change so that’s not a good measure of what is going on.
>>
>> From: Eloi Ribeiro <mail at eloiribeiro.eu>
>> Sent: Tuesday, March 19, 2024 6:56 AM
>> To: PostGIS Users Discussion <postgis-users at lists.osgeo.org>
>> Subject: RE: How to properly rotate a raster
>>
>> Hi all,
>>
>> Here an example with the code.
>>
>> If you compare rast_rotated_from_raster with rast_rotated_from_vector, the latest is rotated as expected and the first is not.
>>
>> I guess, I must be doing something wrong when creating rast_rotated_from_raster.
>>
>> Any guess?
>>
>> -- non rotated raster
>>
>> DROPTABLEIFEXISTS public.rast_not_rotated;
>>
>> CREATETABLEpublic.rast_not_rotated AS
>>
>> SELECT ST_AsRaster(t.geom,1000,1000,'8BUI',1,0) AS rast
>>
>> FROM (SELECT ST_GeomFromText('POLYGON((166989.36462303242 505764.77048559673,
>>
>> 166996.25788353998 505576.476318614,
>>
>> 167265.54134611748 505585.94966528506,
>>
>> 167258.67088940175 505774.2438322678,
>>
>> 166989.36462303242 505764.77048559673))',28992) geom) t;
>>
>> ALTERTABLE public.rast_not_rotated ADD COLUMN rid SERIALPRIMARY KEY;
>>
>> CREATEINDEXrast_not_rotated_rast_gistON public.rast_not_rotated USING gist (st_convexhull(rast));
>>
>> SELECT AddRasterConstraints('public'::name, 'rast_not_rotated'::name,'rast'::name);
>>
>> -- rotate previous raster
>>
>> DROPTABLEIFEXISTS public.rast_rotated_from_raster;
>>
>> CREATETABLEpublic.rast_rotated_from_raster ASSELECT ST_SetRotation(rast, -0.035162353342492736) rast FROM public.rast_not_rotated;
>>
>> ALTERTABLE public.rast_rotated_from_raster ADD COLUMN rid SERIALPRIMARY KEY;
>>
>> CREATEINDEXrast_rotated_from_raster_rast_gistON public.rast_rotated_from_raster USING gist (st_convexhull(rast));
>>
>> SELECT AddRasterConstraints('public'::name, 'rast_rotated_from_raster'::name,'rast'::name);
>>
>> -- rotate polygon and then save as raster
>>
>> DROPTABLEIFEXISTS public.rast_rotated_from_vector;
>>
>> CREATETABLEpublic.rast_rotated_from_vector AS
>>
>> SELECT ST_AsRaster(t.geom,1000,1000,'8BUI',1,0) AS rast
>>
>> FROM (SELECT ST_Rotate(
>>
>> ST_GeomFromText('POLYGON((166989.36462303242 505764.77048559673,
>>
>> 166996.25788353998 505576.476318614,
>>
>> 167265.54134611748 505585.94966528506,
>>
>> 167258.67088940175 505774.2438322678,
>>
>> 166989.36462303242 505764.77048559673))',28992)
>>
>> ,-0.035162353342492736, 166989.36462303242, 505764.77048559673) geom) t;
>>
>> ALTERTABLE public.rast_rotated_from_vector ADD COLUMN rid SERIALPRIMARY KEY;
>>
>> CREATEINDEXrast_rotated_from_vector_rast_gistON public.rast_rotated_from_vector USING gist (st_convexhull(rast));
>>
>> SELECT AddRasterConstraints('public'::name, 'rast_rotated_from_vector'::name,'rast'::name);
>>
>> Eloi
>>
>> On Monday, March 18th, 2024 at 19:11, Eloi Ribeiro <mail at eloiribeiro.eu> wrote:
>>
>>> Hi Regina and all,
>>>
>>> Elaborating a bit more. I have a vector layer consisting of a rectangle (farm plot), and this rectangle is subdivided into several (polygons) stripes parallel to the outer edge of the rectangle. I want to produce a raster where the rectangle and stripes long edges are parallel to the latitude line.
>>>
>>> Using ST_Rotate(geom, -0.035, x, y) in all vector layers together with ST_AsRaster(), I managed to get what I want. Good.
>>>
>>> Nevertheless, if I first produce the raster and then try to rotated it using ST_SetRotation(rast, -0.035), I barely see any rotation. In this step, I must be doing something wrong.
>>>
>>> I would like to first rasterize vector and then rotate the raster. Instead of rasterize non rotated vector and then, again rasterize rotated vector.
>>>
>>> The change in size (width and height), I mentioned in my first message, was while testing with 0.5 rad, not with -0.035 rad.
>>>
>>> -- input raster
>>>
>>> SELECT ST_Width(rast), ST_Height(rast), ST_SkewX(rast) FROM plot_3;
>>>
>>> -- st_width | st_height | st_skewx
>>>
>>> -- ----------+-----------+----------
>>>
>>> -- 1000 | 1000 | 0
>>>
>>> -- output raster of ST_SetRotation(rast, -0.035162353342492736)
>>>
>>> SELECT ST_Width(rast), ST_Height(rast), ST_SkewX(rast) FROM plot_3_rotated_from_raster;
>>>
>>> -- st_width | st_height | st_skewx
>>>
>>> -- ----------+-----------+----------------------
>>>
>>> -- 1000 | 1000 | 0.006952538312005269
>>>
>>> -- output raster of ST_SetRotation(rast, -0.035162353342492736)
>>>
>>> SELECT ST_Width(rast), ST_Height(rast), ST_SkewX(rast) FROM plot_3_rotated_from_vect;
>>>
>>> -- st_width | st_height | st_skewx
>>>
>>> -- ----------+-----------+----------
>>>
>>> -- 977 | 964 | 0
>>>
>>> -- output raster of ST_SetRotation(rast, -0.5)
>>>
>>> SELECT ST_Width(rast), ST_Height(rast), ST_SkewX(rast) FROM plot_3_rotated_from_raster_05;
>>>
>>> -- st_width | st_height | st_skewx
>>>
>>> -- ----------+-----------+---------------------
>>>
>>> -- 1000 | 1000 | 0.09481479675190897
>>>
>>> Cheers,
>>>
>>> Eloi
>>>
>>> On Monday, March 18th, 2024 at 14:01, Regina Obe <lr at pcorp.us> wrote:
>>>
>>>> By size, what do you mean exactly? It would change the width and height.
>>>>
>>>> What are you expecting rotation to do. Perhaps you can give example out of the below and some sample
>>>>
>>>> ST_Width, ST_Height, ST_SkewX before and after and what you were expecting.
>>>>
>>>> https://postgis.net/docs/en/RT_ST_SetRotation.html
>>>>
>>>> From: Eloi Ribeiro <mail at eloiribeiro.eu>
>>>> Sent: Monday, March 18, 2024 4:57 AM
>>>> To: PostGIS Users Discussion <postgis-users at lists.osgeo.org>
>>>> Subject: How to properly rotate a raster
>>>>
>>>> Hi all,
>>>>
>>>> I need to rotate a raster, for that I'm using the function ST_SetRotation, like so:
>>>>
>>>> CREATE TABLE plot_3_rotated AS
>>>> SELECT ST_SetRotation(rast, 0.03516235334249185) rast
>>>> FROM plot_3;
>>>>
>>>> But I see that is not producing the expected results and is changing size of the raster. What am I doing wrong?
>>>>
>>>> Cheers,
>>>>
>>>> Eloi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20240323/e1a4ae96/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 29210 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20240323/e1a4ae96/attachment.png>