6

Is it possible using any of the command line GDL tools such as gdalwarp, etc to perform a rotation of DEM data - a GeoTiff for example?

Anyone know how to do this?


UPDATE1

After clearing out GeoTiff tags and loading my plain Tiff into QGIS with a modified TFW file - I am able to successfully view and export my rotated Tiff DEMs. From what I can tell, the elevation resolution is preserved! However, upon export I decided to view them in MicroDEM and noticed some blocky artifacts. I went back to QGIS and saw the same thing:

MicroDEM enter image description here

QGIS enter image description here

Is this perhaps an error in my transform calculations? I did a -10 degree rotation - Started with this:

0.0000925926
0.0000000000
0.0000000000
-0.0000925926
-157.0956018521
20.9603240742

Ended up with this:

0.00009118591
-0.00001607853
-0.00001607853
-0.00009118591
-157.0956018521
20.9603240742

Math look ok there? What else could this be?

Have you seen this before? Can this not be avoided when performing such a rotation?

(Should I move this updated section to a new gis.stackexchange.com question?)


UPDATE2

I performed the same process with a few additional rotation angles, specifically 45° and 90°. The blocky artifacts seem much less severe to the point that they almost don't appear at all. However, another more obvious problem arises at these rotations. The data appears to brighten and become more white overall at certain rotations. In other words, the elevation rises generally in the entire model and detail is therefore lost. First, notice the original view in QGIS:

Original (no rotation), 45° and 90° - respectively:

Original 45 degress 90 degrees

In fact, the same problem is occurring with my original -10 degree rotation (though more subtle):

-10 degrees

This leads me to believe that if the data is rotated at all then I will get this clipping effect. Remember, this is the data as QGIS presents it upon opening the .TIFF file. I have not exported anything or applied any special processes beyond modifying the .TFW file.

What in the world could be going on here? How can I avoid these distortions when rotating the data?


UPDATE3

As noted in the comments, I resolved the first issue using the bi-linear re-sampling. The second issue was resolved by setting the min/max correctly.

  • 1
    It depends totally on the software that is taking care of rendering or other processing of the "rotated" DEM. If you fiddle with .tfw you do not touch a single pixel in the tiff file, you can even make it read-only. If the result is not good you can make a new permanently rotated image with gdalwarp gdal.org/gdalwarp.html gdalwarp rotated_with_tfw.tif warped physically.tif Test different resampling options but nearest or bilinear should be good for DEM. Read also about the other gdalwarp options. – user30184 Jan 25 '15 at 9:45
  • @user30184 Either my math is wrong, or I am missing some fundamental step to properly rotating my data. Any tips you can provide are appreciated. – dtmland Jan 26 '15 at 4:52
  • What do you mean with the clipping effect? The top left pixel of the original image is nailed to a fixed coordinate with any rotation angle. The original rectangular image needs a bigger rectangular box to fit in when it is rotated, despite at 90, 180, and 270 degree, and the corner areas are filled with black no-data pixels. Would you rather see the corners filled with white pixels? – user30184 Jan 26 '15 at 6:22
  • The clipping effect is the "higher elevation" portions of the terrain becoming all white and apparently losing most of its detail? It is most obvious in the 45 degree rotation shown above. Do you see that occurrence? – dtmland Jan 26 '15 at 6:30
  • 2
    You have to compare the rendering of the grayscale. QGIS sets limits to 2% and 98% of the min and max values by default. That might get changed when you rotate and introduce NODATA values. You can set the min and max values for rendering manually to compare it better. – AndreJ Jan 26 '15 at 6:54
3

You can create a local CRS with an oblique mercator projection, and transform the data with gdalwarp and gdal_translate into it. See my advice here:

Using customized Coordinate System in ArcGIS Desktop?

This should work with 16-bit or grayscale data the same way. Paletted colours shoud be expanded to RGBA in advance.


UPDATE

Using QGIS, create custom CRS without (tmerc) and with (omerc) the rotation, using the same coordinates for the center of projection:

+proj=tmerc +lat_0=51.4 +lon_0=7 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
+proj=omerc +lat_0=51.4 +lonc=7 +alpha=-10 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

In a first run, use Raster -> Projections -> Warp on your tif to the rotation point:

gdalwarp -overwrite -s_srs EPSG:4326 -t_srs "+proj=tmerc +lat_0=51.4 +lon_0=7 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -of GTiff D:/Download/N51E007.hgt D:/Download/51_7.tif

In QGIS, you will not see a difference, because the pixel values have not changed. Then use Raster -> Conversion -> Translate:

gdal_translate -a_srs "+proj=omerc +lat_0=51.4 +lonc=7 +alpha=-10 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -of GTiff D:/Download/51-7.tif D:/Download/51-7rot.tif

and the result looks like expected:

enter image description here

Make sure to copy the style from the unrotated to the rotated layer to see the same styling.

  • Your answer in the linked question - it references a solution using the arcgic-desktop software? Should the steps be relatively the same for QGIS? (I realize I didn't specify this in my original question) – dtmland Jan 25 '15 at 8:01
  • The question was for arcgis, but my answer was done in QGIS. – AndreJ Jan 25 '15 at 9:08
  • Sorry I was mixed up - I'm mostly new to this stuff. I'm going over your post - and struggling to derive how exactly the degree of rotation is specified? Is it implied when providing the appropriate coordinates? – dtmland Jan 26 '15 at 1:45
  • 2
    The rotation is given with +alpha=, clockwise in degrees. – AndreJ Jan 26 '15 at 5:17
  • See my extended answer for clarification. – AndreJ Jan 26 '15 at 7:13
2

Rather simple method is to write a new world file (.tfw) which contains rotation parameters. You can make such with OpenOffice Calc, for example. If you have a GeoTIFF file which contains reoreferencing info as stored into the image tag you must clear the geotiff tags and create a baseline tiff to start with. It can be done with gdal_translate:

gdal_translate -of GTiff -co profile=baseline -co tfw=yes geotiff.tif baseline.tif

Example of an original, not rotated tfw. Zeroes as the value of (B) and (C) tell that there is no rotation.

0.0092605453   (A)
0.0000000000   (B)
0.0000000000   (C)
-0.0092605453  (D)
34.0044913837  (E)
38.9955086163  (F)

Use the following formulas for calculating new values into rotated .tfw file. If rotation angle "r" is given in degrees, then (angle)=r*Pi/180

(A')=A*Cos(angle)
(B')=A*Sin(angle)
(C')=-D*Sin(angle)
(D')=D*Cos(angle)
(E')=E
(F')=F 

Test with rotating a sample tif with -10 degrees (negative value means clockwise rotation.

New .tfw file will be

0.0091198568
-0.0016080768
-0.0016080768
-0.0091198568
34.0044913837
38.9955086163

With your tiff file (e.g. baseline.tif) and rotated .tfw (e.g. baseline.tfw) in the same directory - open the tiff file into QGIS. You should see an image rotated around the top-left pixel.

enter image description here

Internally the rotated .tfw file converts into GDAL GeoTransform array. The same result without fiddling with .tfw file can be achieved by setting GeoTransform programmatically for example from a Python script.

gdalinfo rotated.tif
Driver: GTiff/GeoTIFF
Files: rotated.tif
       rotated.tfw
       rotated.tif.aux.xml
Size is 216, 216
Coordinate System is `'
GeoTransform =
  34.0007354937, 0.009119856799999999, -0.0016080768
  39.0008725831, -0.0016080768, -0.009119856799999999
  • Thanks for the details - I believe I've just about got it but have run into a snag. I've updated the question. – dtmland Jan 25 '15 at 7:48

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.