Writing CRS to Raster

[Philip Hestvik]

Hi John,

I am working on a simple script to create a raster of an height above ground model from lidar data. My lidar is in ~1km square tiles and is in a projected coordinate system. 

wbe = wbw.WbEnvironment()

wbe.verbose = True # Let each of the function calls output to stdout.

print(wbe.version()) # Let's see what version of WbW we're working with

fp_in = Path(r'D:\343_5934.laz')

fp_out = Path(rf'D:\{fp_in.stem}_HAG.tif')

if not fp_out.parent.exists():

    fp_out.parent.mkdir(parents=True, exist_ok=True)

lidar =  wbe.read_lidar(str(fp_in))

wkt = lidar.get_well_known_text()

crs = CRS.from_user_input(wkt)

lidar_hag = wbe.height_above_ground(input=lidar)

dem_hag = wbe.lidar_tin_gridding(input_lidar=lidar_hag, interpolation_parameter = "elevation", returns_included = "all", cell_size = 0.5)

dem_hag = wbe.set_nodata_value(raster=dem_hag, back_value=-9999.0)

dem_hag = wbe.fill_missing_data(dem_hag, filter_size=11) # Notice that we overwrite `dem` here.

wbe.write_raster(dem, str(fp_out), compress=True) # Compression is good, but it is a bit slower so here we won't use it.

If I do this though, the raster is written and the values are as expected, but it is missing a CRS.

I looked into the dem_hag.configs.epsg_code to see if I could "set it" but I didn't see that option.

I saw that I could create a copy the config to a variable and set attributes this way:

out_configs = dem_hag.configs

out_configs.epsg_code = epsg_code(crs)[0] #calls the function and returns the EPSG code

out_configs.nodata = -9999.0

dem = wbe.new_raster(out_configs)

From here I can create a new raster with the CRS, but I then need to copy the raster values over.

Is there an existing function I am not seeing to solve this problem?

Thanks,
Philip

[Philip Hestvik]

Here is my solution using the rioxarray library

        wbe.write_raster(dem_hag, str(fp_out), compress=True)

        try:

            with rioxarray.open_rasterio(fp_out) as xds:

                xds.rio.set_crs(f"EPSG:{epsg_code(crs)[0]}")

                #xds.rio.crs

                fp_out_crs = fp_out.parent.joinpath(f"{fp_out.stem}_crs.tif")

                xds.rio.to_raster(fp_out_crs)

            fp_out.unlink()

            fp_out_crs.rename(fp_out)

        except:

            print(f'{fp_out} failed to write CRS')

If there is a way to write the CRS with WBW I am all ears, I just couldn't figure it out myself

-Philip