Aller au contenu

Data To Raster

Command

Bases: BaseCommand

Source code in back/iarbre_data/management/commands/data_to_raster.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
class Command(BaseCommand):
    help = "Convert Data polygons to raster."

    def add_arguments(self, parser):
        """Add arguments to the command."""
        parser.add_argument(
            "--grid-size", type=int, default=5, help="Grid size in meters"
        )

    def handle(self, *args, **options):
        output_dir = str(BASE_DIR) + "/media/rasters/"
        resolution = 1
        grid_size = options["grid_size"]
        kernel_size = grid_size / resolution
        if not kernel_size.is_integer():
            raise ValueError(
                f"Grid size should be a multiple of the resolution: {resolution}m"
            )
        kernel_size = int(kernel_size)
        # remove test data
        all_cities_union = City.objects.exclude(code="38250").aggregate(
            union=Union("geometry")
        )["union"]
        minx, miny, maxx, maxy = all_cities_union.extent

        width = int((maxx - minx) / resolution)
        height = int((maxy - miny) / resolution)
        transform = from_origin(minx, maxy, resolution, resolution)

        width_out = int((maxx - minx) / grid_size)
        height_out = int((maxy - miny) / grid_size)
        transform_out = from_origin(minx, maxy, grid_size, grid_size)

        for factor_name in FACTORS.keys():
            log_progress(f"Processing {factor_name}")
            if factor_name == "QPV":
                continue
            rasterize_data_across_all_cities(
                factor_name,
                height,
                width,
                height_out,
                width_out,
                transform,
                transform_out,
                all_cities_union=all_cities_union,
                grid_size=kernel_size,
                output_dir=output_dir,
            )

add_arguments(parser)

Add arguments to the command.

Source code in back/iarbre_data/management/commands/data_to_raster.py 
95
96
97
98
99
def add_arguments(self, parser):
    """Add arguments to the command."""
    parser.add_argument(
        "--grid-size", type=int, default=5, help="Grid size in meters"
    )

rasterize_data_across_all_cities(factor_name, height, width, height_out, width_out, transform, transform_out, all_cities_union, grid_size=5, output_dir=None)

Convert Data polygons to a single binary raster across all cities to avoid border effects.

This function rasterizes the geometries of a specified factor across all cities, applies a convolution to aggregate the raster values into larger blocks, and saves the resulting raster to a file.

Parameters:

Name Type Description Default
factor_name str

Name of the factor to transform to raster.

 required
height int

Height of the output raster.

 required
width int

Width of the output raster.

 required
height_out int

Height of the output raster after convolution.

 required
width_out int

Width of the output raster after convolution.

 required
transform Affine

Affine transformation for the factor transformation.

 required
transform_out Affine

Affine transformation for the raster output.

 required
all_cities_union GEOSGeometry

GEOSGeometry containing the union of all city geometries.

 required
grid_size int

Size of the convolution kernel. Defaults to 5.

 5
output_dir str

Directory to save the raster file. Defaults to None.

 None

Returns:

Type Description
None

None
Source code in back/iarbre_data/management/commands/data_to_raster.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def rasterize_data_across_all_cities(
    factor_name: str,
    height: int,
    width: int,
    height_out: int,
    width_out: int,
    transform: rasterio.Affine,
    transform_out: rasterio.Affine,
    all_cities_union: GEOSGeometry,
    grid_size: int = 5,
    output_dir: str = None,
) -> None:
    """
    Convert Data polygons to a single binary raster across all cities to avoid border effects.

    This function rasterizes the geometries of a specified factor across all cities, applies a convolution
    to aggregate the raster values into larger blocks, and saves the resulting raster to a file.

    Args:
        factor_name (str): Name of the factor to transform to raster.
        height (int): Height of the output raster.
        width (int): Width of the output raster.
        height_out (int): Height of the output raster after convolution.
        width_out (int): Width of the output raster after convolution.
        transform (rasterio.Affine): Affine transformation for the factor transformation.
        transform_out (rasterio.Affine): Affine transformation for the raster output.
        all_cities_union (GEOSGeometry): GEOSGeometry containing the union of all city geometries.
        grid_size (int, optional): Size of the convolution kernel. Defaults to 5.
        output_dir (str, optional): Directory to save the raster file. Defaults to None.

    Returns:
        None
    """
    max_count = grid_size * grid_size
    os.makedirs(output_dir, exist_ok=True)

    qs = Data.objects.filter(factor=factor_name, geometry__intersects=all_cities_union)
    factor_df = load_geodataframe_from_db(qs, [])
    log_progress(f"Rasterizing {factor_name}")
    # Rasterize the shapes
    raster = rasterize(
        factor_df.geometry,
        out_shape=(height, width),
        transform=transform,
        fill=0,
        default_value=1,
        dtype=np.uint8,
    )
    if len(raster[raster > 0]) == 0:
        raise ValueError(f"{factor_name} is producing a blank tif.")
    log_progress("Sum on 5x5")
    kernel = np.ones((grid_size, grid_size))
    coarse_raster = ndimage.convolve(raster, kernel, mode="constant", cval=0)[
        0 : height_out * grid_size : grid_size,
        0 : width_out * grid_size : grid_size,
    ]
    coarse_raster = (coarse_raster / max_count * 100).astype(np.int8)
    # Save the raster to file
    output_path = os.path.join(output_dir, f"{factor_name}.tif")
    log_progress(f"Saving {factor_name}")
    with rasterio.open(
        output_path,
        "w",
        driver="GTiff",
        height=height_out,
        width=width_out,
        count=1,
        dtype=raster.dtype,
        crs="EPSG:2154",
        transform=transform_out,
    ) as dst:
        dst.write(coarse_raster, 1)