Obtain coordinates and corresponding pixel values from GeoTiff using python gdal and save them as numpy array

Question

How can I obtain projected coordinates as well as the actual pixel values at those coordinates from a GeoTiff file and then save them into a numpy array? I have arsenci020l.tif file, and its coordinates are in meters. Below is the abridged output of gdalinfo I ran on it.

~$ gdalinfo arsenci020l.tif 
Driver: GTiff/GeoTIFF
Files: arsenci020l.tif
 arsenci020l.tfw
Size is 10366, 7273
Coordinate System is:
PROJCS["Lambert Azimuthal Equal Area projection with arbitrary plane grid; projection center 100.0 degrees W, 45.0 degrees N",
 GEOGCS["WGS 84",
 DATUM["WGS_1984",
 SPHEROID["WGS 84",6378137,298.257223563,
 AUTHORITY["EPSG","7030"]],
 AUTHORITY["EPSG","6326"]],
 PRIMEM["Greenwich",0],
 UNIT["degree",0.0174532925199433],
 AUTHORITY["EPSG","4326"]],
 PROJECTION["Lambert_Azimuthal_Equal_Area"],
 PARAMETER["latitude_of_center",45],
 PARAMETER["longitude_of_center",-100],
 PARAMETER["false_easting",0],
 PARAMETER["false_northing",0],
 UNIT["metre",1,
 AUTHORITY["EPSG","9001"]]]
Origin = (-6086629.000000000000000,4488761.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
...

There was a similar question here about obtaining lat/long coordinates from tiff (Obtain Latitude and Longitude from a GeoTIFF File) and the answer showed how to obtain only top left x and y pixel coordinates. I need to obtain ALL projected pixel coordinates as well as get the pixel values and save them in a numpy array. How can I do it?

• Mike T: You want 10366 × 7273 = over 75 million points?
• irakhman: @MikeT I think so,I don't really know of a better solution of how to approach the problem I'm trying to solve:I need to find the closest pixel coordinate from this dataset to each centroid of US block and then assign the corresponding pixel value to that block.From searching around I realized that cKDTree query is going to help me with nearest neighbor search.Python function for the algorithm asks for a "tree" to query as numpy array.In order to make a "tree" out of all pixel coordinates from this dataset,I need to store all of them somehow.If you have a better solution,please let me know!
• geotheory: I cannot understand why this procedure requires so much code in Python.

Answer 1

This should get you going. The raster values are read using rasterio, and pixel centre coordinates are converted to Eastings/Northings using affine, which are then converted to Latitude/Longitude using pyproj. Most arrays have the same shape as the input raster.

import rasterio
import numpy as np
from affine import Affine
from pyproj import Proj, transform
fname = '/path/to/your/raster.tif'
# Read raster
with rasterio.open(fname) as r:
 T0 = r.transform # upper-left pixel corner affine transform
 p1 = Proj(r.crs)
 A = r.read() # pixel values
# All rows and columns
cols, rows = np.meshgrid(np.arange(A.shape[2]), np.arange(A.shape[1]))
# Get affine transform for pixel centres
T1 = T0 * Affine.translation(0.5, 0.5)
# Function to convert pixel row/column index (from 0) to easting/northing at centre
rc2en = lambda r, c: (c, r) * T1
# All eastings and northings (there is probably a faster way to do this)
eastings, northings = np.vectorize(rc2en, otypes=[float, float])(rows, cols)
# Project all longitudes, latitudes
p2 = Proj(proj='latlong',datum='WGS84')
longs, lats = transform(p1, p2, eastings, northings)

• mitchus: When using this, I get the message "AttributeError: 'DatasetReader' object has no attribute 'affine'" for the line "T0 = r.affine"
• Translunar: @mitchus Apparently affine is just an alias for transform, and the alias has been removed from the most recent version of rasterio. I edited the answer but it looks like it needs to be peer-reviewed since I'm new here. :)
• Translunar: It also looks like the indexes are wrong for A.shape, which has only two dimensions.
• bw4sz: rasterio pads a 1 channel dimension so that shape is always (channel, row, col)

Answer 2

would add as comment, but a bit long - in case you wanted to use gdal/ogr within python - something like this might work (hacked together from some other code i had - not tested!) This also assumes that rather than finding the nearest raster pixel to a polygon centroid, you simply query the raster at the xy of the centroid. i have no idea what the speed tradeoff might be...

from osgeo import gdal,ogr
fc='PathtoYourVector'
rast='pathToYourRaster'
def GetCentroidValue(fc,rast):
 #open vector layer
 drv=ogr.GetDriverByName('ESRI Shapefile') #assuming shapefile?
 ds=drv.Open(fc,True) #open for editing
 lyr=ds.GetLayer(0)
 #open raster layer
 src_ds=gdal.Open(rast) 
 gt=src_ds.GetGeoTransform()
 rb=src_ds.GetRasterBand(1)
 gdal.UseExceptions() #so it doesn't print to screen everytime point is outside grid
 for feat in lyr:
 geom=feat.GetGeometryRef()
 mx=geom.Centroid().GetX()
 my=geom.Centroid().GetY()
 px = int((mx - gt[0]) / gt[1]) #x pixel
 py = int((my - gt[3]) / gt[5]) #y pixel
 try: #in case raster isnt full extent
 structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_Float32) #Assumes 32 bit int- 'float'
 intval = struct.unpack('f' , structval) #assume float
 val=intval[0]
 except:
 val=-9999 #or some value to indicate a fail
 feat.SetField('YOURFIELD',val)
 lyr.SetFeature(feat)
 src_ds=None
 ds=None
GetCentroidValue(fc,rast)