Generating input data

[Peter9192]

This issue is to track work on generating input data. Doesn't exactly fit in WRF code, but it is relevant for using the functionality devised in #4

[Peter9192]

I tried producing input data using convert_geotiff, using the following command:

convert_geotiff -m -999 -s 0.01 -d "Building albedo" HERITAGE-albedo_raster.tif

However, the -999 used to signal to convert_geotiff does not seem to get picked up as expected. The input tiff has -999 as NODATA value, but cells without data are converted as all others (divided by scale factor --> -99900) and stored as int16 --> to large, so we get rubbish values.

Need to find an alternative way to write the data correctly.

[Peter9192]

I now figured out how to convert the geotiff to WPS input format using plain python:

import numpy as np
import rasterio
import matplotlib.pyplot as plt

SCALE_FACTOR = 0.01
MISSING_VALUE = -999.0
TILE_BDR = 3
TILE_SIZE = 100

# endianness = 'little' if np.little_endian else 'big'


# Read tiff file
with rasterio.open("HERITAGE-albedo_raster.tif") as src:
    source = src.read(1)
    nodata = src.nodata
    crs = src.crs
    nx, ny = src.shape


# Build output array (including halo)
shape = (nx + 2 * TILE_BDR, ny + 2 * TILE_BDR)
dtype = "<i2"  # little endian | signed integer | 2 bytes
dest = np.ones(shape) * MISSING_VALUE


# Fill with data
# Apply scale factor and deal with missing values
scaled_source = np.where(source == nodata, MISSING_VALUE, source / SCALE_FACTOR)
dest[TILE_BDR:-TILE_BDR, TILE_BDR:-TILE_BDR] = np.flipud(scaled_source)

filename = f"{1:05d}-{nx:05d}.{1:05d}-{ny:05d}"
dest.astype(dtype).tofile(filename)

This doesn't tile the data as convert_geotiff does, however, that is also not strictly necessary. I can develop this further to also set the CRS from the geotiff file, and automatically generate the index file with correct endianness etc.

To verify this worked, I tried to be as close as possible to what convert_geotiff does (including e.g. flipud and a tile boundary, which is not strictly necessary I think). This leads to the following verification:

tile_shape = (106, 106)
ref = np.fromfile("00001-00100.00001-00100", dtype="<i2").reshape(tile_shape)
new = np.fromfile("00001-00225.00001-00493", dtype="<i2").reshape(shape)[:106, :106]

ref_fixed_missing = np.where(ref < 0, -999, ref)

fig, axs = plt.subplots(1, 3)
axs[0].imshow(ref_fixed_missing)
axs[1].imshow(new)
axs[2].imshow(new - ref_fixed_missing)

axs[0].set_title("convert_geotiff")
axs[1].set_title("custom python")
axs[2].set_title("diff")

fig.savefig("verify.png", bbox_inches="tight")

The only difference is that the reference from convert_geotiff has some values inside the boundary, which is not supposed to be the case if I understand it correctly.

[Peter9192]

This version is quite decent and has been tested through geogrid:

import numpy as np
import rasterio

# Configuration (boring)
MISSING_VALUE = -999
SCALE_FACTOR = 0.01
dtype = "<i2"  # little endian, signed integer, wordsize 2

# Read GeoTIFF
with rasterio.open("HERITAGE-albedo_raster.tif") as src:
    data = src.read(1).astype(np.float32)
    nodata = src.nodata
    ny, nx = data.shape
    transform = src.transform

    # Grid spacing (degrees)
    dx = transform.a
    dy = -transform.e

    # Lower-left grid cell center
    lon, lat = src.xy(0, 0)


# Scale to integer space
scaled_data = np.where(data == nodata, MISSING_VALUE, data / SCALE_FACTOR)

# Write binary
filename = f"{1:05d}-{nx:05d}.{1:05d}-{ny:05d}"
scaled_data.astype(dtype).tofile(filename)

# Write index
index_file_content = f"""projection = regular_ll
known_x    = 1
known_y    = {ny}
known_lat  = {lat:.6f}
known_lon  = {lon:.6f}
dx         = {dx:.5e}
dy         = {dy:.5e}
tile_x     = {nx}
tile_y     = {ny}
tile_z     = 1
type = continuous
signed = yes
wordsize = 2
endian = little
row_order = top_bottom
tile_bdr = 0
missing_value = {MISSING_VALUE}
scale_factor = {SCALE_FACTOR}
units = "NO UNITS"
description = "Building albedo"
"""

with open("index", "w") as f:
    f.write(index_file_content)


# check output: plot the binary file
import matplotlib.pyplot as plt

wps_binary = np.fromfile(filename, dtype="<i2").reshape((ny, nx))
wps_masked = np.ma.masked_values(wps_binary, MISSING_VALUE)
fig, ax = plt.subplots()
im = ax.imshow(wps_masked, cmap="viridis")
plt.colorbar(im, orientation="horizontal")
fig.savefig("wps_binary.png", bbox_inches="tight")