Yo can access to the centroid of each pixel as doc says with ST_PixelAsCentroids (postgis 2.1)
SELECT x, y, val, ST_AsText(geom) FROM (SELECT (ST_PixelAsCentroids(rast, 1)).* FROM dummy_rast WHERE rid = 2) foo;
x | y | val | st_astext
---+---+-----+--------------------------------
1 | 1 | 253 | POINT(3427927.775 5793243.975)
2 | 1 | 254 | POINT(3427927.825 5793243.975)
Now you have the geom , and its trivial to split the lat/lon into separate columns if you want.
To get the appropriate cell coordinates from your latitude and longitude you need to know the coverage of the netCDF file (for example, the coordinate of the upper left cell and the size of the cells in the x and y directions). @Luke 's solution below works for a srtaight x/y raster with no cell rotation. If you don't have that, you could look at the affine
library to transform latitude and longitude to cell coordinates. Given GDAL GeoTransform in the form:
- Upper left cell
x
coordinate
x
cell size in projected units
x
directional rotation (usually 0)
- Upper left cell
y
coordinate
y
directional rotation (usually 0)
- Negative
y
cell size in projected units
So for example the -237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0
would translate as an upper left cell of the raster with coordinates -237481.5, 237536.4
, and a 425.0
unit square with no rotation.
Using the affine
library you can transform this into a transformation object like so:
from affine import Affine
aff = Affine.from_gdal(-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)
Which can transform cell coordinates to projected coordinates:
xs = np.array([0, 1, 2])
ys = np.array([3, 4, 5])
x_coords, y_coords = aff * (xs, ys)
Or (in your case) from coordinates back to cell coordinates:
xs, ys = ~aff * (np.array(lons), np.array(lats))
These values are floats, so you'll need to transform them into integers to get cell coordinates you can use.
xs = np.round(xs).astype(np.int)
ys = np.round(ys).astype(np.int)
You can then use these as indexes to your netCDF4 array (use the latest version of netCDF4 - 1.2.1 as this means you don't need to sort or remove duplicates).
variable = fnc.variables['variable'][xs, ys] # be careful of your index order
This will return a square array due to the slightly different netCDF indexing, but you can get the actual values that you're after as the diagnoal:
values = variable.diagonal()
Best Answer
GeoServer will serve the raster value as a GetFeatureInfo of the WMS service. To see this in action, go to the Layer-Preview page, select your layer and just click on the map. On the bottom, you will see your raster value (either RGB triplet or your TIFF channel(s)).
Below is an example from http://projects.bryanmcbride.com/leaflet/wms_info.html . This link expire but the internet archive yields the javascript code below
There is also an OpenLayers example.