1

My goal is to get the coordinates of all pixels with a specified value (ex: maximum and minimum pixel values).

# imports
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio

# load tif
src = rasterio.open('./data/thhz.tif')

# read band
band = src.read(1)

# get maximum and minimal values
max = np.max(band)
min = np.min(band)

# convert raster into shapefile
# got this from https://gis.stackexchange.com/questions/346288/extract-all-pixels-values-from-geotiff-with-python
### this process is too slow from here need to use other methods
px_vals = []

for x in range(band.shape[0]):
    for y in range(band.shape[1]):
        px_vals.append({'x': x, 
                        'y': y,
                        'value': band[x, y]})

# convert list to dataframe
vals = pd.DataFrame(data=px_vals, columns=['x', 'y', 'value'])

# convert dataframe to geopandas 
vals['geometry'] = gpd.points_from_xy(vals['x'], vals['y'])
vals = gpd.GeoDataFrame(vals, geometry='geometry', crs='EPSG:4326')

# transform coordinates with affine of raster
# get affine matrix with src.profile
vals['geometry'] = vals.affine_transform([0.001, 0.0, 2.616837667, 0.0, -0.001, 11.815923996])
# this also was not well converted

# slice values with coordinates
min_pixels = vals[vals['value'] == max]
max_pixels = vals[vals['value'] == min]

2 Answers 2

1

Looping over a numpy array will always be slower than using vectorized numpy operations only. Try:

min_rows, min_cols = np.where(band == np.min(band))
max_rows, max_cols = np.where(band == np.max(band))

band == np.min(band) gives back a boolean array, broadcasted to be in the same shape as band, and containing True where the element is equal to the min value.

np.where() gives back indices where the condition is true, if you only specify the first argument.

5
  • Is it possible to get the coordinates instead of #row & #cols?
    – pyaj
    Mar 26, 2022 at 0:22
  • min_lons, min_lats = transform.xy(thhz.transform, min_rows, min_cols)
    – pyaj
    Mar 26, 2022 at 11:39
  • max_lons, max_lats = transform.xy(thhz.transform, max_rows, max_cols)
    – pyaj
    Mar 26, 2022 at 11:40
  • from rasterio import transform
    – pyaj
    Mar 26, 2022 at 11:40
  • Please add these to your solution. Thanks
    – pyaj
    Mar 26, 2022 at 11:41
0

I just worked out a solution

def coordinates_and_values(raster, pixel_values):
    """
    this function extracts longitude and latitudes of list of raster pixel values
    raster: load raster with rasterio
    pixel_values: list of pixel values 
    exception: pixel value must be found more than one time in raster else len of float error
    """
    df = pd.DataFrame()                         # create dataframe
    raster_band = raster.read(1)                # read raster band
    for i in pixel_values:                      # iterate between list of pixel values
        rows, cols = np.where(raster_band == i) # extract row and column numbers for each pixel
        rows, cols = transform.xy(raster.transform, rows, cols) # transform row and column numbers to coordinates
        values = np.array([i] * len(rows))      # create array containing n pixel value of n coordinates
        df_i = pd.DataFrame(zip(rows, cols, values), columns=['lon','lat', 'pixel']) # create dataframe for one pixel value
        df = df.append(df_i)                    # append to get dataframe of lon and lat of list of pixel values
    return df

Example:

coordinates_and_values(src, [-7, -6, -5, 6])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.