Which points in projected raster corresponds to certain latitude,longitude location?

Question

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:

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?

Answer 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)

Answer 2

If you have a source image with these limits

φ_min, λ_min and φ_max, λ_max

that has been warped to target image with these limits

N_min, E_min and N_max, E_max

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

Just rearrange

(φ - φ_min) / (φ_max - φ_min) = (N - N_min) / (N_max - N_min)

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.