It is not clear precisely what form the data are in, but ultimately any solution will have to use the equations for the Mercator projection of an ellipsoid. It sounds like it's possible to identify the latitude and longitude for at least two pixels, such as a given one and one at an origin. Along with the datum, this will suffice.
The datum, among other things, describes the size and shape of the ellipsoid. We will need its eccentricity e, which is related to the flattening via
e^2 = f * (2 - f).
For WGS 84, for instance,
1 / f = 298.257 223 563
giving e^2 = 0.00669438 and e is approximately 0.08181919084. (For the sphere, f and e are both zero, considerably simplifying the following equations.)
We also need the equatorial radius (semi-major axis) of the ellipsoid, a; for WGS 84, a = 6,378,137.0 meters.
Every image or map has a (global) scale: it is the amount by which the projection coordinates are uniformly multiplied to place projected points (which are in meters) on the images (which is in pixels). Let's call this scale S; its unit of measurement is therefore pixels per meter. Then, relative to an origin on the image, the x (easting) and y (northing) coordinates of the Mercator projection of the point with longitude lambda and latitude phi are
x = S * a * lambda
y = S * a * [log(tan(pi/4 + phi/2)) + (e/2) * (log(1 - e * sin(phi)) - log(1 + e * sin(phi)))]
Although this may look nasty, notice that e and a are known constants, S is an unknown constant, and phi is known for each pixel (presumably, it is the latitude of the pixel's center). Let's encapsulate that information by abbreviating the formula for y as
y = S * a * M(phi; e, a)
(Incidentally, M(phi; 0, a) = log(tan(pi/4 + phi/2)) has a particularly simple form for the sphere.)
Suppose we have a pixel at row i0 and column j0 corresponding to (lambda0, phi0) on the earth, and another pixel at row i1 and column j1 corresponds to (lambda1, phi1). Then we can deduce S in two ways from the two equations. The x equation will work provided the pixels are in different columns and the y equation will work provided the pixels are in different rows:
S = (j1 - j0) / (lambda1 - lambda0) / a
S = (i1 - i0) / (M(phi1; e, a) - M(phi0; e, a)) / a
At least one of these equations will work. With S now known, we can compute the Mercator coordinates for any point (lambda, phi) on the globe directly. We can also invert the projection numerically or by using a (rapidly) iterative algorithm or with a trigonometric series, as described by Snyder on p. 44.
From the information in the question we might guess that S equals 1 pixel / 1 Km = 0.001 pixels per meter. However, that presumes the projection has true scale at the Equator. The allusion to a reference latitude in a comment to the question suggests that the projection has already been rescaled by some amount. That is why it appears necessary to compute S. Otherwise, we would not need information about two pixels: the latitude, longitude, row, and column of a single pixel would suffice.
Reference
Snyder, John P. Map Projections--A Working Manual. USGS professional paper 1395, 1987.
You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module.
An example:
from math import floor
from osgeo import gdal,ogr
import struct
src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'
src_ds=gdal.Open(src_filename)
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
rb=src_ds.GetRasterBand(1)
ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY() #coord in map units
#Convert from map to pixel coordinates.
px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
px = floor(px) #x pixel
py = floor(py) #y pixel
structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)
print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
Alternatively, since the reason you gave for not using numpy was to avoid reading the entire array in using ReadAsArray(), below is an example that uses numpy and does not read the entire raster in. It uses the built-in gdal.ApplyGeoTransform() function in order to deal with axes rotations.
from math import floor
from osgeo import gdal,ogr
src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'
src_ds=gdal.Open(src_filename)
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
rb=src_ds.GetRasterBand(1)
ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY() #coord in map units
#Convert from map to pixel coordinates.
px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
px = floor(px) #x pixel
py = floor(py) #y pixel
intval=rb.ReadAsArray(px,py,1,1)
print intval[0] #intval is a numpy array, length=1 as we only asked for 1 pixel value
Best Answer
To find which feature the coordinate is near, you first need to build an R-Tree index of bounding boxes or envelope of each feature. A popular library for this is libspatialindex.
Secondly, you would then need to know for each of the matched features from your R-tree, which ones match. GDAL/OGR does have some operations based on GEOS to see if the geometry "contains" the point of interest, and you can extract the field info. See the OGRGeometry Class Reference for
Contains,Touches,Within, etc. that can perform the geometry relation operator.Personally I wouldn't do any of the above, because it would take me too long to figure out the coding, and I know there are faster ways. (I would loading all shapefiles in a PostGIS database, then query the locations using SQL.) However, it is understandable if this GIS functionality needs to be embedded in an existing project without adding more complicated dependencies, such as a database server.