carbon.py

from pathlib import PurePath
import subprocess as sp
from parallelpipe import stage
import logging
import glob
import sys


@stage(workers=96)
def calc(cabon_tiles, area_tiles):
    for carbon in cabon_tiles:

        tile_id = "_".join(PurePath(carbon).parts[-1].split("_")[:2])

        area = None
        for a in area_tiles:
            if tile_id in a:
                area = a
                break

        if not area:
            logging.warning("Could not find matching area tile for: " + tile_id)
            break

        output = "output/{}_mt_co2_pixel_2000".format(tile_id)

        cmd = ["gdal_calc.py", "-A", carbon, "-B", area, "--outfile", output, "--calc","A*B/10000/2*44/12",
               "--NoDataValue", "0","--type", "Float32", "--co", "NBITS=16", "--co", "TILED=YES",
               "--co", "COMPRESS=LZW", "--co", "PREDICTOR=3", "--co", "SPARSE_OK=TRUE"]

        try:
            logging.info("Processing tile: " + tile_id)
            sp.check_call(cmd)
        except sp.CalledProcessError:
            logging.warning("Failed to process tile: " + tile_id)
        else:
            yield output


def main():

    root = logging.getLogger()
    root.setLevel(logging.DEBUG)

    handler = logging.StreamHandler(sys.stdout)
    formatter = logging.Formatter(
        "%(asctime)s - %(name)s - %(levelname)s - %(message)s"
    )
    handler.setFormatter(formatter)

    root.addHandler(handler)

    carbon_tiles = (glob.glob("carbon/*.tif"))
    area_tiles = (glob.glob("area/hanson_2013*.tif"))

    pipe = carbon_tiles | calc(area_tiles)

    for output in pipe.results():
        logging.info("Processed : " + output)

    logging.info("Done")


if __name__ == "__main__":
    main()

gadm.py

import subprocess as sp
import logging
from parallelpipe import stage


def get_longitude(x):
    if x >= 0:
        return str(x).zfill(3) + "E"
    else:
        return str(-x).zfill(3) + "W"


def get_latitude(y):
    if y >= 0:
        return str(y).zfill(2) + "N"
    else:
        return str(-y).zfill(2) + "S"


def bbox():
    for lat in range(-90,90,10):
        for lon in range(-180,180,10):
            xmin = lon
            xmax = lon + 10

            ymin = lat
            ymax = lat + 10

            yield xmin, ymin, xmax, ymax


@stage(workers=46)
def rasterize(tiles):

    for tile in tiles:
        xmin = tile[0]
        ymin = tile[1]
        xmax = tile[2]
        ymax = tile[3]

        top = get_latitude(ymax)
        left = get_longitude(xmin)
        tile_id = "{}_{}".format(top, left)

        output = "gadm36_uid_{}.tif".format(tile_id)

        cmd = ["gdal_rasterize", "-a", "UID", "-te", str(xmin), str(ymin), str(xmax), str(ymax),
               "-tr", "0.00025", "0.00025", "-tap", "-ot", "UInt32", "--co", "NBIT=19",
               "--co", "TILED=YES", "--co", "COMPRESS=LZW", "--co", "PREDICTOR=2",
               "--co", "SPARSE_OK=TRUE", 'PG:"user=postgres dbname=postgres"', "gadm", output]

        try:
            sp.check_call(cmd)
        except sp.CalledProcessError:
            logging.warning("Failed to process tile: " + tile_id)
        else:
            yield output


def main():

    root = logging.getLogger()
    root.setLevel(logging.DEBUG)

    handler = logging.StreamHandler(sys.stdout)
    formatter = logging.Formatter(
        "%(asctime)s - %(name)s - %(levelname)s - %(message)s"
    )
    handler.setFormatter(formatter)

    root.addHandler(handler)

    tiles = [tile for tile in bbox()]

    pipe = tiles | rasterize()

    for output in pipe.results():
        logging.info("Processed : " + output)

    logging.info("Done")


if __name__ == "__main__":
    main()