[GIS] Calculate pixel coordinates for each pixel in raster? EPSG:27700 (British National Grid)

coordinate systempixelpythonraster

I have some Met Office rain radar files (NIMROD binary files). These contain raster data where each pixel is a 1km grid square. I am trying to calculate the geographic position of each of these pixels so that I have the coordinates of each position recorded in a separate file. Using this I can then programmatically work out which catchment (polygon) the pixel belongs in.

Broadly I am asking how best to do this.

I have some information about the data. It is a grid 2175 by 1725. Each pixel is a 1km square. I am given the bounding box which is in EPSG:27700 (British National Grid) projection. The values are apparently -404500 1549500 1320500 -625500. These seemed strange to me as they appear to be outside the BNG bounds?

start northing 1549500.0, start easting -404500.0

Is it too simplistic to simply take the distance between the two easting bounds (for example) and decide by the number of pixels in the row? This should locate the centre of the most far west pixel if I am not mistaken. Then it would be a simple matter to move through the row. I want to end up with a matrix of the coordinates for each pixel.

Alternatively would it be possible to take the most east bound and somehow add 1km west (the result being given as the new easting value)?

The ultimate idea is that I have a catchment polygon within the UK and I want to identify each raster pixel which falls over this area or otherwise intersects the area. I have several of these polygons. I don't want to use a GIS tool as I don't need to visualise it. I just need to know what pixels intersect which polygons.

You can see above my initial thoughts about how to go about this but if my approach is too simplistic due to the nature of project please let me know.

Best Answer

The raster bounds, size and resolution you state are not consistent and I think you may need to reinvestigate your data. For example, if it really is a 1km pixel size then your raster spans 2,175km; that is rather large and well beyond recommended usage for BNG. Also your given bounding box is similarly not consistent with a 1km resolution.

Generally speaking however, from a rectangular raster one may infer pixel coordinates as follows. In this case I assume zero-based indexing with respect to the NW corner:

Given pixel array index i,j:

E = E0 + i * s + s/2
N = Nm - j * s - s/2

Where:

  • E0 is minimum Easting;
  • Nm is maximum Northing;
  • s is pixel size/resolution.

Note that the s/2 terms ensure that the derived coordinate is the centre of the cell.

Related Question