From “Raster-vector interactions”

Last updated: 2023-02-25 13:39:06

---

Exercise 11-c

import numpy as np
import geopandas as gpd
import rasterio
import rasterio.mask
from rasterio.plot import show
# Prepare raster
src = rasterio.open('output/carmel.tif')
# Prepare polygons
stat = gpd.read_file('data/statisticalareas_demography2019.gdb')
stat = stat[['STAT11', 'SHEM_YISHUV', 'geometry']]
sel = stat['SHEM_YISHUV'] == 'חיפה'
stat = stat[sel]
stat = stat.to_crs(src.crs)
# Mask
out_image, out_transform = rasterio.mask.mask(src, stat['geometry'].to_list(), crop=True, nodata=-9999)
out_image = out_image[0]
# To 'float'
out_image1 = out_image.astype(np.float64)
out_image1[out_image1 == -9999] = np.nan

Exercise 11-e

import numpy as np
import pandas as pd
import shapely.geometry
import geopandas as gpd
import rasterio
import rasterio.mask
from rasterio.plot import show
# Prepare raster
src = rasterio.open('output/carmel.tif')
r = src.read(1)
# Prepare polygons
stat = gpd.read_file('data/statisticalareas_demography2019.gdb')
stat = stat[['STAT11', 'SHEM_YISHUV', 'geometry']]
sel = stat['SHEM_YISHUV'] == 'חיפה'
stat = stat[sel]
stat = stat.to_crs(src.crs)
# Raster to polygons
shapes = rasterio.features.shapes(r, mask=r > 525, transform=src.transform)
pol = list(shapes)
geom = [shapely.geometry.shape(i[0]) for i in pol]
geom = gpd.GeoSeries(geom, crs=src.crs)
values = [i[1] for i in pol]
values = pd.DataFrame({'elev': values})
result = gpd.GeoDataFrame(data=values, geometry=geom)
# Dissolve
result = result.unary_union

glue("exercise-11-e", result)

result.area

421200.0

Exercise 11-d

import numpy as np
import geopandas as gpd
import rasterio
import rasterio.mask
from rasterio.plot import show
# Prepare raster
src = rasterio.open('output/carmel.tif')
# Prepare polygons
stat = gpd.read_file('data/statisticalareas_demography2019.gdb')
stat = stat[['STAT11', 'SHEM_YISHUV', 'geometry']]
sel = stat['SHEM_YISHUV'] == 'חיפה'
stat = stat[sel]
stat = stat.to_crs(src.crs)
# Mask
out_image, out_transform = rasterio.mask.mask(src, stat['geometry'].to_list(), crop=False, nodata=-9999)
out_image = out_image[0]
# To 'float'
out_image1 = out_image.astype(np.float64)
out_image1[out_image1 == -9999] = np.nan

Exercise 11-f

import numpy as np
import rasterio
import shapely.geometry
import rasterstats
import geopandas as gpd

# Prepare raster
src = rasterio.open('output/carmel.tif')
r = src.read(1)
r1 = r.copy()
r1[r > 525] = np.arange(0, r[r > 525].size)

# Raster to polygons
shapes = rasterio.features.shapes(r1, mask=r > 525, transform=src.transform)
pol = list(shapes)
geom = [shapely.geometry.shape(i[0]) for i in pol]
geom = gpd.GeoSeries(geom, crs=src.crs)

# Polygons to points
geom = geom.centroid
result = gpd.GeoDataFrame(geometry=geom)

# Extract values to points
x = rasterstats.point_query(
    result, 
    r, 
    nodata = src.nodata, 
    affine = src.transform,
    interpolate='nearest'
)
result['elev'] = x

# Plot
result.plot(column='elev', legend=True);

Exercise 11-g

import pandas as pd
import geopandas as gpd
import rasterio
from rasterstats import zonal_stats

# Prepare raster
src = rasterio.open("output/carmel.tif")
r = src.read(1)
r[(r != -9999) & (r < 100)] = 0
r[r >= 100] = 1

# Prepare polygons
stat = gpd.read_file('data/statisticalareas_demography2019.gdb')
stat = stat[['SHEM_YISHUV', 'STAT11', 'geometry']]
sel = stat['SHEM_YISHUV'] == 'חיפה'
stat = stat[sel]
stat = stat.reset_index(drop = True)
stat = stat.to_crs(src.crs)

# Extract
result = zonal_stats(stat, r, nodata=-9999, affine=src.transform, stats='mean')
result = pd.DataFrame(result)
stat1 = pd.concat([stat, result], axis=1)

# Plot
stat1.plot(column="mean", legend=True);

Exercise 11-h

import pandas as pd
import geopandas as gpd
import rasterio
from rasterstats import zonal_stats

# Prepare raster
src = rasterio.open('output/sentinel2.tif')
r = src.read()
r = r / 10000

# Prepare polygons
pol = gpd.read_file('data/muni_il.shp')
pol = pol[['Muni_Eng', 'geometry']]
sel = ["Be'er Sheva", "Lehavim", "Omer", "Lakiye", "Segev Shalom"]
pol = pol[pol["Muni_Eng"].isin(sel)]
pol = pol.to_crs(src.crs)

# Extract
result = []
for band in range(r.shape[0]):
    tmp = zonal_stats(pol, r[band], nodata=-9999, affine=src.transform, stats='mean')
    tmp = pd.DataFrame(tmp)
    tmp = tmp.rename(columns = {'mean': 'band_' + str(band)})
    tmp.index = sel
    result.append(tmp)
result = pd.concat(result, axis=1)

# Plot
result.T.plot();