2

What I have is a SRTM 3 arc-second raster file in WGS84 geographic coordinate system (here is that file). I then project it to the Azimuthal Equidistant projection with GDAL:

gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_ts=40.81266 +lon_0=14.414252" -r near -of GTiff C:/vesuvius2_wgs84.tif C:/vesuvius2_aeqd.tif

I would like to know, how can I determine which point in the newly projected raster (C:/vesuvius2_aeqd.tif) corresponds to certain latitude,longitude location?

For example I would like to check where does the location of 14.4, 40.8 lay in the "C:/vesuvius2_aeqd.tif" file:

enter image description here

And for example it lays at -1051, -450 projected coordinates.

How can I calculate this?

EDIT:

Is it possible to know the exact pixel which corresponds to the required geographic coordinate location? For example, the 14.4, 40.8 geographic coordinate, corresponds to the pixel: 200x210 in the projected raster. Is that possible?

  • Why? Warping is destructive, it sounds like you just need point value extraction - do you really need to remodel this raster? – mdsumner Apr 30 '16 at 12:24
  • @mdsumner, please forgive me if I explained it in a bad way. Yes, I just need a point transformation, not transformation of the whole raster. Is it possible to provide the exact pixel which corresponds to the given latitude,longitude location? For example: it is the pixel 200x210? – marco Apr 30 '16 at 12:52
  • Are you asking What GDAL/OGR function does point-to-point map projection? – Martin F Apr 30 '16 at 19:30
  • Also, "+proj=merc" != "+proj=aeqd" – mkennedy May 2 '16 at 18:12
  • @MartinF, I believe yes. mkennedy: I did not mention Mercator projection in my reply. What is the meaning of your reply? – marco May 4 '16 at 20:59
1

For an example of how to do this in code (Python/GDAL), look at this example in the GDAL cookbook, in particular the world2Pixel() function. I've added some comments to show what the variables are:-

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0] # minimum x geospatial
  ulY = geoMatrix[3] # minimum y geospatial
  xDist = geoMatrix[1] # width of cell in geospatial
  yDist = geoMatrix[5] # height of cell in geospatial coords
  rtnX = geoMatrix[2] # rotation (not used)
  rtnY = geoMatrix[4] # rotation (not used)
  pixel = int((x - ulX) / xDist) # x pixel coord
  line = int((ulY - y) / xDist) # y pixel coord (Bug?)
  return (pixel, line)

Actually, looking at that, I'm puzzled by why xDist is used in the function to assign to line, rather than yDist.. I suspect this is a bug which might bite if cells aren't square.

There are some edge cases where this needs to be more complicated, such as rasters with rotation where north is not 'up' (these are shown as rtnX and rtny in the code but never used)

It's worth having a look through the various other code examples on that page, as they also include examples of the inverse operation (getting geospatial coordinates from pixel coordinates).

EDIT

I've found there's actually a GDAL example script called val_at_coord.py (here's the source). If you call it with something like

python -m val_at_coord -display_xy -coordtype=georef long lat filename.tiff

.. it should return the coords.

On my system at least this is in a zip file and not installed by default. You might need to hunt it down extract this somewhere and run it from a terminal window (or OSGeo4W shell on Windows)

  • thank you Steven. Is it possible to do this, without using the programming (or python)? Is it possible to do it directly through GDAL command prompt? – marco May 10 '16 at 10:26
  • I've found an example gdal script which seems to do the job, have edited the answer. – Steven Kay May 10 '16 at 20:14
1

If you have a source image with these limits

φmin, λmin and φmax, λmax

that has been warped to target image with these limits

Nmin, Emin and Nmax, Emax

surely you can do simple bilinear interpolation to map a general point from φ, λ to N, E?

Just rearrange

(φ - φmin) / (φmax - φmin) = (N - Nmin) / (Nmax - Nmin)

and similarly for λ and E.

However, if you wish to use GDAL/OGR for point-to-point map projection, consult the section Coordinate Transformation of the gdal.org/osr_tutorial:

The OGRCoordinateTransformation class is used for translating positions between different coordinate systems. New transformation objects are created using OGRCreateCoordinateTransformation(), and then the OGRCoordinateTransformation::Transform() method can be used to convert points between coordinate systems.

  • thank you for the reply. Is it possible to do this, without using the programming? Is it possible to do it directly through GDAL command prompt? – marco May 10 '16 at 10:26

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.