0

I'm currently using this code to convert a raster file to a geodataframe:

import rasterio 
from rasterio.features import shapes
mask = None 

with rasterio.open(#INSERT TIF FILE HERE) as src:
    image = src.read(1) # first band, not sure yet how to do it with multiple bands
    results = (
    {'properties': {'raster_val': v}, 'geometry': s}
    for i, (s, v) 
    in enumerate(
        shapes(image, mask=mask))) geoms = list(results)

import geopandas as gpd
gpd_polygonized_raster = gpd.GeoDataFrame.from_features(geoms)

The problem is, the geodataframe is showing upside down instead of its expected oreintation.

Any help on this would be appreciated. Thank you!

Take note that the TIFF file has a projection already of EPSG:4326.

Oshean Lee Garonita
  • 33
  • 1
  • 5

2 Answers2

0
  • you can use shapely affine_transform()
  • have picked up a sample GEOTIFF* to make this working example
import rasterio
from rasterio.features import shapes
import geopandas as gpd
from shapely.affinity import affine_transform as T
from pathlib import Path
import plotly.express as px
import requests
from pathlib import Path

# https://www.sciencebase.gov/catalog/item/53f5a87ae4b09d12e0e8547b
url = "https://www.sciencebase.gov/catalog/file/get/53f5a87ae4b09d12e0e8547b?f=__disk__7f%2Fe7%2F94%2F7fe7943c77f2c6c4eb4c129153fd4e80f6079091"
f = Path.cwd().joinpath("FAA_UTM18N_NAD83.tif")
r = requests.get(url, stream=True, headers={"User-Agent": "XY"})
with open(f, "wb") as fd:
    for chunk in r.iter_content(chunk_size=128):
        fd.write(chunk)

results = []

with rasterio.open(f) as src:
    crs = src.crs
    for layer in range(1,20):
        try:
            image = src.read(layer)  # use all bands / layers
            results += [
                {"properties": {"raster_val": v, "layer":layer}, "geometry": s}
                for i, (s, v) in enumerate(shapes(image, mask=src.dataset_mask()))
            ]
        except IndexError as e:
            print(e)
            break


gdf = gpd.GeoDataFrame.from_features(results, crs=crs)

# flip top and bottom using affine transform
gdf["geometry"] = gdf["geometry"].apply(lambda g: T(g, [1, 0, 0, -1, 0, 0]))

# exclude areas where it's white in source TIFF
gdf.to_crs("EPSG:4326").loc[gdf["raster_val"].lt(255)].plot(column="raster_val")

enter image description here

Rob Raymond
  • 29,118
  • 3
  • 14
  • 30
  • Hi Rob! Thanks for this. I will try! But just for context, I'm currently using this image: https://drive.google.com/file/d/1bcaYoVltRxswWiOuEL_QhqrNag-8FkrX/view?usp=sharing The file is already at epsg:4326. – Oshean Lee Garonita Dec 02 '21 at 17:11
  • For affine transformation, I am not sure on what to put after T – Oshean Lee Garonita Dec 02 '21 at 17:14
  • it's documented here: https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations I tried with your tif 1) I don't believe it is EPSG:4326 as co-ordinates are not valid for EPSG:4326. 2) to plot needed to modify last line to: `gdf.loc[gdf["raster_val"].gt(0)].plot(aspect=1, column="raster_val")` – Rob Raymond Dec 02 '21 at 18:59
  • Hmm, that's weird. When I opened the file on a GIS software (QGIS), it's showing the right location. Anyway I'll try doing affine transformation. – Oshean Lee Garonita Dec 03 '21 at 15:13
0
you could vectorize the tif file using gdal_polygonize and then read the vector file with geopandas
import subprocess
import geopandas as gpd

# execute 'gdal_polygonize.py' in the shell
cmd = ['gdal_polygonize.py', tif_path, vector_path]
subprocess.call(cmd, shell=True)

# read the vector file (.shp, .geojson, etc..)
gpd_polygonized_raster =  gpd.read_file(vector_path)
EOEAGLE
  • 1