Use Solweig without Qgis dependencies

[paolodep36]

Are there any ways to use Solweig model to run analysis on a python environment such as colab? I don't want to rely on qgis. Let's imagine i want to calculate the mean radiant temperature and the Pet and i have all the data to do that, is there any workaround to just achieve this?

[j3r3m1]

I think that @b-kode has previously done something like that.

[paolodep36]

@biglimp yes DSM_KRbig.tif is still in 3007. If i write the same code but using rasterio, coordinates are set well `import rasterio
from rasterio.warp import transform_bounds, transform
import numpy as np

def process_dsm_with_rasterio(dsm_path, feedback=None):
with rasterio.open(dsm_path) as src:
# Read the DSM data
dsm = src.read(1).astype(float)

    # Handle no-data values
    nd = src.nodata
    dsm[dsm == nd] = 0.0

    # Adjust for negative values
    if dsm.min() < 0:
        dsmraise = np.abs(dsm.min())
        dsm += dsmraise
        if feedback:
            feedback(f'Digital Surface Model (DSM) included negative values. DSM raised with {dsmraise}m.')

    # Extract CRS and bounds
    old_crs = src.crs
    bounds = src.bounds

    # Print the CRS info
    print("Old CRS:", old_crs)

    # Define the WGS 84 coordinate system
    new_crs = 'EPSG:4326'

    # Transform the bounds to the new CRS
    transformed_bounds = transform_bounds(old_crs, new_crs, *bounds)
    
    print("Transformed Bounds (WGS 84):", transformed_bounds)

    # Calculate the center point of the bounds
    center_x = (bounds.left + bounds.right) / 2
    center_y = (bounds.bottom + bounds.top) / 2

    # Transform the center point
    center_lon, center_lat = transform(old_crs, new_crs, [center_x], [center_y])
    
    print("Center (WGS 84):", center_lon[0], center_lat[0])

    # Calculate scale
    scale = 1 / src.transform[0]

    # Calculate altitude
    alt = np.median(dsm)
    if alt < 0:
        alt = 3

    # Feedback
    if feedback:
        feedback(f'Longitude derived from DSM: {center_lon[0]}')
        feedback(f'Latitude derived from DSM: {center_lat[0]}')

    return center_lon[0], center_lat[0], alt, scale

Example usage

def simple_feedback(message):
print(message)

dsm_path = '/content/UMEP/Data _goteborg/Goteborg_SWEREF99_1200/DSM_KRbig.tif'
longitude, latitude, altitude, scale = process_dsm_with_rasterio(dsm_path, feedback=simple_feedback)
print(f'Longitude: {longitude}, Latitude: {latitude}, Altitude: {altitude}, Scale: {scale}')`

[nilswallenberg]

Yes, I had similar issues, and they were related to gdal version. If you try to derive latitude and longitude as in gdal 2.x but with gdal 3.x or vice versa you will get wrong coordinates. Maybe this is the reason?

[paolodep36]

I was working on colab just to make the process as generic as possible. Changing gdal version on a machine turns out to be an often not immediate and very delicate operation. I set lat and lon according to the txt produced by solweig on qgis version and imposing a gdal version = 2. the results are very similar, the difference is a few tenths of a celsius degree. The gdal version of qgis I am using is 3.2.1 while in colab it is 3.6.4. In any case the code is very sensitive to latitude and longitude. It would be interesting to understand where lat and lon are used and if they are recalculated in some other process related to solweig

[paolodep36]

I came up with a solution using this code but i don't know if i am doing wrong with the remaining logic of the code ` gdal_dsm = gdal.Open(filepath_dsm)

    # Handle no-data values
    dsm = gdal_dsm.ReadAsArray().astype(float)
    nd = gdal_dsm.GetRasterBand(1).GetNoDataValue()
    dsm[dsm == nd] = 0.
    if dsm.min() < 0:
        dsmraise = np.abs(dsm.min())
        dsm = dsm + dsmraise
        print(f'Digital Surface Model (DSM) included negative values. DSM raised with {dsmraise}m.')
    else:
        dsmraise = 0

    # Get the original CRS
    dsm_ref = gdal_dsm.GetProjection()
    old_cs = osr.SpatialReference()
    old_cs.ImportFromWkt(dsm_ref)
    print("Original CRS (EPSG):", old_cs.GetAttrValue("AUTHORITY", 1))

    # Define the WGS 84 coordinate system
    new_cs = osr.SpatialReference()
    new_cs.ImportFromEPSG(4326)
    print("New CRS (EPSG):", new_cs.GetAttrValue("AUTHORITY", 1))

    # Reproject the DSM to EPSG:4326 using gdal.Warp
    warp_options = gdal.WarpOptions(dstSRS='EPSG:4326', format='MEM')
    reprojected_dsm = gdal.Warp('', gdal_dsm, options=warp_options)

    # Get geotransform and dimensions of the reprojected DSM
    reprojected_geotransform = reprojected_dsm.GetGeoTransform()
    minx = reprojected_geotransform[0]
    maxx = minx + reprojected_geotransform[1] * reprojected_dsm.RasterXSize
    miny = reprojected_geotransform[3] + reprojected_geotransform[5] * reprojected_dsm.RasterYSize
    maxy = reprojected_geotransform[3]

    # Calculate the center point of the bounds
    lon = (minx + maxx) / 2
    lat = (miny + maxy) / 2

    # Calculate the median altitude
    alt = np.median(dsm)
    if alt < 0:
        alt = 3

    print(f'Longitude derived from DSM: {lon}')
    print(f'Latitude derived from DSM: {lat}')

`

[paolodep36]

It seems that i had success to solve the little differences between qgis result and the generic code. I noticed that the test dataset of goteborg have no data value not properly set. I set it all to -9999 and now results are the same. Also the cdsm is not good georeferenced!.Hope this could help someone else in the future!