From “Raster-vector interactions”#

Last updated: 2023-05-27 16:40:07


Exercise 11-c#

import numpy as np
import geopandas as gpd
import rasterio
from rasterio.plot import show
# Prepare raster
src = rasterio.open('output/carmel.tif')
# Prepare polygons
stat = stat[['STAT11', 'SHEM_YISHUV', 'geometry']]
sel = stat['SHEM_YISHUV'] == 'חיפה'
stat = stat[sel]
stat = stat.to_crs(src.crs)
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
from rasterio.plot import show
# Prepare raster
src = rasterio.open('output/carmel.tif')
# Prepare polygons
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
from rasterio.plot import show
# Prepare raster
src = rasterio.open('output/carmel.tif')
# Prepare polygons
stat = stat[['STAT11', 'SHEM_YISHUV', 'geometry']]
sel = stat['SHEM_YISHUV'] == 'חיפה'
stat = stat[sel]
stat = stat.to_crs(src.crs)
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')
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[(r != -9999) & (r < 100)] = 0
r[r >= 100] = 1

# Prepare polygons
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 = r / 10000

# Prepare polygons
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();