7

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.

Improve this question
9
  • 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. Commented Jan 25, 2015 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. Commented Jan 26, 2015 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? Commented Jan 26, 2015 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? Commented Jan 26, 2015 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. Commented Jan 26, 2015 at 6:54

2 Answers 2

Reset to default
4

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.

Improve this answer
5
  • 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) Commented Jan 25, 2015 at 8:01
  • The question was for arcgis, but my answer was done in QGIS. Commented Jan 25, 2015 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? Commented Jan 26, 2015 at 1:45
  • 2
    The rotation is given with +alpha=, clockwise in degrees. Commented Jan 26, 2015 at 5:17
  • See my extended answer for clarification. Commented Jan 26, 2015 at 7:13
3

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
Improve this answer
1
  • Thanks for the details - I believe I've just about got it but have run into a snag. I've updated the question. Commented Jan 25, 2015 at 7:48

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.