Finding the centroid of every pixel in a raster Python

Question

I am trying to get the centroid of every pixel in a aviris hyperspectral image (a raster) but my results are not quite as expected.

The centroids don't seem to line up with the center of every pixel, as shown by the images below. In fact, it's quite smaller than the aviris image.

Here is the code I used to get the centroids. Based off of this related answer

aviris_path = 'aviris/f190802t01p00r18_rfl_v1l1/f190802t01p00r18_corr_v1l1_img'
def get_centroid_raster(img_path):
read_img = rasterio.open(img_path)

#read in the image to get the shape
open_img = read_img.read().transpose(1,2,0)

#number of rows and columns
num_rows = open_img.shape[0]
num_cols = open_img.shape[1]

#hold the resulting lon and lat
hold_centroid_coordinates = np.zeros((num_rows, num_cols, 2))

for row in range(num_rows):
    for col in range(num_cols):

        the_coords = rasterio.transform.xy(
            read_img.transform,
            row,
            col,
            offset = 'center'
        )

        hold_centroid_coordinates[row, col, 0] = the_coords[0]
        hold_centroid_coordinates[row, col, 1] = the_coords[1]

#flatten the lon and lat into a 1-dimensional array
longitude = hold_centroid_coordinates[:,:,0].flatten()
latitude = hold_centroid_coordinates[:,:,1].flatten()

#put the longitude and latitude into a dataframe
coordinate_dataframe = pd.DataFrame({'longitude': longitude, 'latitude' : latitude})
return coordinate_dataframe

Not sure where I'm going wrong.

Answer 1

You can accomplish this without the nested for loops by using np.meshgrid(). I think this will speed up your processing time as well, give this script a try:

import rasterio
import numpy as np
import geopandas as gpd
from shapely.geometry import Point
path = 'your_raster.tif'
with rasterio.open(path) as src:
    band1 = src.read(1)
    height = band1.shape[0]
    width = band1.shape[1]
    cols, rows = np.meshgrid(np.arange(width), np.arange(height))
    xs, ys = rasterio.transform.xy(src.transform, rows, cols)
    lons = np.array(xs)
    lats = np.array(ys)
points = gpd.GeoSeries(
    list(zip(lons.flatten(), lats.flatten()))).map(Point)

# use the feature loop in case shp is multipolygon
geoms = points.values
features = [i for i in range(len(geoms))]

out = gpd.GeoDataFrame(
    {'feature': features, 'geometry': geoms}, crs=src.crs)
out.to_file("pixel_center_points.shp")

input:

output:

Credit to these previous posts that helped me come up with the solution:

Get coordinates of all pixels in a raster with rasterio

How to write Shapely geometries to shapefiles?

Making shapefile from Pandas dataframe?