ThinkHazard - EarthQuake Processing

eq_global.py

import logging
import sys
import numpy
import rasterio
import fiona
from rasterio import features
from shapely.geometry import shape
from shapely.ops import transform
from functools import partial
import pyproj

logging.basicConfig(stream=sys.stderr, level=logging.INFO)
logger = logging.getLogger('rasterize_geometry')

threshold = 98.655

return_periods = ['0250', '0475', '1000']
rasters = []



# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.drivers():
    for rp in return_periods:
        logger.info('Reading data for return period %s' % rp)

        with rasterio.open('eq%sfpgm.tif' % rp) as src:
            src_data = src.read()
            profile = src.profile

        dst_data = (src_data > threshold).astype(rasterio.uint8)

        profile.update(
            dtype=rasterio.uint8,
            count=1,
            nodata=0)
        rasters.append(dst_data)

    logger.info('Merge data (total)')
    total = numpy.zeros(rasters[0].shape)
    for rp in rasters:
        total += rp

logger.info('Reading admin divisions shapefile')
collection = fiona.open("administrativedivision.shp")

logger.info('Rasterizing each feature of the admin divisions')
for f in list(collection):
    code = f['properties']['CODE']
    project = partial(
            pyproj.transform,
            pyproj.Proj(init='epsg:3857'),
            pyproj.Proj(init='epsg:4326')
            )
    geom = transform(project, shape(f['geometry']))
    division = features.rasterize(
            [(geom, 1)],
            out_shape=src.shape,
            transform=src.transform)

    masked = numpy.ma.masked_array(total, mask=~division.astype(bool))

    max_ = numpy.max(masked)

    logger.info('Level for %s is %s'
            % (f['properties']['NAME'], numpy.max(masked)))