I have a lot of radar images that come with image coordinates (824x648). The projection should be Austria Lambert/ETRS89 (EPSG: 3416).

So to perform the georeferencing I need to convert ROW/COL to LAMBDA/PHI, right? I got 6 parameters that are apparently used for this purpose (v0 thru v5: 8.194 50.437 0 0 1018.18 1018.18), but how am I doing this? As I have around 1000 radar images, this should be done somehow automated, so I guess this is a affine transformation issue, right (e.g. 6-tuple geotransform in GDAL)? How do I figure out what is what for my 6 parameter transformation without having any additional informtion? There are different types of parameter order around apparently.

  • Are those 6 numbers the coordinates of a corner in lat-long and the cell heights and widths in metres in that coordinate system? Is Austria 8W and 50N and about 824km by 648km? Any chance of (a link to) a sample image? – Spacedman Sep 14 '16 at 13:23
  • All other information I can provide is the projection origin (Long/Lat: 13.333333° 47.0°) at pixel position x/y: 358 387; standard parallels: 46.0° 49.0°; corner coordinates upper left pixel: 8.194° 50.437°. Unfortunately I can't provide a sample image due to legal reasons, sorry :-( – GeoEki Sep 14 '16 at 13:29

Assuming your 6 numbers are top left bounds in EPSG:4326 coordinates, and the last two figures are the cell width and height in meters in the EPSG:3416 system then the procedure is as follows.

  1. Find the top-left corner in EPSG:3416 coordinates.
  2. Add width * (pixel rows+1) to get the eastern limit
  3. Subtract height * (pixel columns + 1) to get the southern limit
  4. We now know the EPSG:3416 coordinates of the grid system.
  5. Use gdal_translate to create an intermediate VRT with control points mapping pixels to EPSG:3416 coords
  6. Use gdalwarp to do the transformation.

So you need the gdal command line tools - you didn't mention any GIS tools so I used them.

I've used a test image which I scaled to 824x648 pixels, called rose.png. The commands were then:

 gdal_translate -of VRT -a_srs EPSG:3416 \
   -gcp 0 0 34787.89 738735.2 \
   -gcp 824 648 874786.4 78954.57 \
   -gcp 0 648 34787.89 78954.57 \
   rose.png rose3416c.vrt

 gdalwarp -s_srs EPSG:3416 -t_srs EPSG:3416 rose3416c.vrt rose3416c.tif

This produces a GeoTIFF that seems to sit nicely over Austria, and looks rectangular on an EPSG:3416 base map:

enter image description here

I used R's projection functions to work out the coordinates of the corners, and a bit of maths. Note this may be half a cell out if the coordinates are cell centres and not edges. Or even more than that out if I've messed up. No warranty. Also, I'm surprised this isn't easier.

If you have 1000 of these to do I'd write a command line script to do it.

Here's how to work out the coordinates using Qgis Python code at the Qgis Python Console:

First set up two projections and the transform object:

>>> atproj = QgsCoordinateReferenceSystem(3416, QgsCoordinateReferenceSystem.PostgisCrsId)
>>> llproj = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.PostgisCrsId)
>>> xform = QgsCoordinateTransform(llproj, atproj)

Now you can convert x-y coordinates from lat-long to Austrian:

>>> xform.transform(QgsPoint(8.194,50.437))

So I can store this point as origin, add the offsets, and get the far corner:

>>> w = 1018.18
>>> h = 1018.18
>>> origin = QgsPoint(xform.transform(QgsPoint(8.194,50.437)))
>>> origin
>>> corner = QgsPoint(origin.x() + 824*w, origin.y() - 648*h)
>>> corner

I'm not sure why the X coordinate is slightly off from the one I computed a similar way in R - possibly a slightly different set of projection parameters or I missed something. Bit its less than 20m off in 1km grid squares so hopefully no big deal.

  • Your solutions does the job amazingly well. Too bad I have no idea how you derived the correct corner coordinates. I just tried to do the job with the affine package, but my values are entirely different than yours. Maybe I'll give it another shot tomorrow, but your solutions works and solves my problem, thank you. – GeoEki Sep 14 '16 at 16:30
  • Added python code for coordinate transformation. – Spacedman Sep 14 '16 at 17:10

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.