14

I'm just starting in the field of GIS (mostly worked with just python before) and got stuck on a problem. I'm looking for a way to polygonize an already masked GeoTIFF raster with population data in a way that I get no connected pixels. This is needed for further processing where each pixel needs to be the same area. I got as far as to extract all shapes and put them as shapes into a GeoDataFrame, as follows:

with rasterio.open('./Paris.tiff') as raster:
    image = raster.read(1)
    crs = raster.crs
    list_pop = [
        {'cell_pop': value, 'geometry': shape(shp)}
        for i, (shp, value) 
        in enumerate(shapes(image, connectivity=0, transform=raster.transform))
        if value > raster.nodata
    ]
    
df = gpd.GeoDataFrame(list_pop, crs=crs).to_crs(epsg=4326)
df.plot()

This does however connect pixels with the same value (connectivity value of rasterio.features.shapes is minimal 4) and result in the following image.

geodataframe with connected pixels

Is there a way to get all pixels out without them being connected like in the center of the image?

3
  • I figured out a temporary (ugly) solution, to add a random() to every pixel so it deviates by a tiny amount and doesn't get co-joined into an extended shape. I get a dirty feeling though, there must be a better way..
    – Ivotje50
    Jun 21, 2020 at 18:04
  • 2
    I just discovered that I have the same problem. I just want the square grid of the original data with no weird polygons from fusing neighboring grids. Clearly this SHOULD be an option of rasterio, but sadly is not. Also sadly, nobody seems to know how to do this because nobody answered your/our question. Oct 8, 2020 at 9:19
  • 1
    Probably interesting question here where goal not to polygonize instead to convert geotiff into pandas dataframe or geopandas dataframe. After the conversion, the polygonisation could be done. It shows a demo with some geotiff format with georasters but cannot get it working with some geotiff files, perhaps the geotiff standard is varying?
    – hhh
    Oct 29, 2020 at 10:43

2 Answers 2

4

exactly this answer:

You can also do it by xarray or rioxarray library.

import rioxarray as rxr
import xarray as xr

dataarray = rxr.open_rasterio('file.tif')
dataarray = xr.open_rasterio('file.tif')

df = dataarray[0].to_pandas()
3

To elaborate on the answer from @Fee, you can use rioxarray and geopandas to achieve your desired output by following these general steps:

  1. Read raster (my example uses multiband tiff as input and outputs first band to .shp)
  2. Compute pixel centroids and convert to geodataframe while preserving intensity values (consider renaming as population in your case)
  3. Convert geometries to geodataframe (and export)

optional: reduce size of buffer so pixels don't touch by setting 'multiplier' parameter < 1

  1. Output geometries to polygons shapefile

Code:

import geopandas as gpd
import numpy as np
import rasterio
import rioxarray as rx

infile = '25cm_raster.tif'

#set this parameter < 1 to reduce size of pixels, i.e., 0.98 for 2% reduction (optional)
multiplier = 1

#function to retrieve buffer resolution from raster metadata
def get_res(filename):
    image = rasterio.open(filename)
    #get resolution as 1/2 average of x and y pixel dimensions
    gt = image.transform
    pixelSizeX = gt[0]
    pixelSizeY =abs(gt[4])
    res = (pixelSizeX + pixelSizeY)/4
    return res
 
if __name__ == "__main__":

    print("Reading raster...")
    raster = rx.open_rasterio(infile)
    print(f'Read successfully. The dimensions of the raster are: {str(raster[0].shape)}.')

    buff = get_res(infile)
    print(buff)

    #get first band of raster
    band = raster[0]
    x, y, intensity = band.x.values, band.y.values, band.values
    x, y = np.meshgrid(x, y)
    x, y, intensity = x.flatten(), y.flatten(), intensity.flatten()

    print("Converting to GeoDataFrame...")

    #create new geoseries with centroid geometries
    centroids = gpd.GeoDataFrame(geometry=gpd.GeoSeries.from_xy(x, y, crs=band.rio.crs))
    centroids['intensity'] = intensity
    
    #create new geoseries with pixel geometries
    pixels = centroids.buffer(buff*multiplier, cap_style=3)
    polygons = gpd.GeoDataFrame(geometry=pixels, crs=band.rio.crs)
    polygons['intensity'] = intensity

    # Saving GeoDataFrames to shapefile
    #centroids.to_file(f"""vectorized_centroids_{infile.split('.')[0]}.shp""", crs=band.rio.crs)
    polygons.to_file(f"""vectorized_pixels_{infile.split('.')[0]}.shp""", crs=band.rio.crs)

Original image:

enter image description here

Centroids:

enter image description here

Output (Polygons):

enter image description here

Optional (Reduced margins):

enter image description here

Credits to this post by spatial dev guru and this answer by franck theeten.

2
  • 1
    Please, don't use shapefiles.... Use GeoPackage... Shapefiles are horribly outdated with quirks such as field name lenght limits... Also GeoPackage can store multiple geometries in one file, while shapefile needs multiple files for one geometry....
    – Fee
    May 31, 2023 at 19:48
  • 1
    Thanks for the insight. Shapefiles have been a quirky mess since I've used them so I can understand the need for a better format.
    – Kartograaf
    Jun 5, 2023 at 20:24

Your Answer

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

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