From 16affa52b4a1c10a21afb20ec134378e22047418 Mon Sep 17 00:00:00 2001 From: KatKatKateryna <89912278+KatKatKateryna@users.noreply.github.com> Date: Mon, 1 Jul 2024 22:06:32 +0800 Subject: [PATCH] 2.19 release_to_MAIN (#196) 2.19 updates: add GisFeature class speed up raster conversions simplify large polygon display Meshes on Send specify QGIS max version generate installers for Mac (we had it already, just needed a change in CI) * add Gis feature class on send and receive (#195) * send GisFeature * remove display val from child geometries * receive GisFeature * fix * ensure units are strings * fix failing test * qgis max version (#194) * raster conversion speed (#193) * draft based on c# (adjust offset/rotation; XY directions) * speed up, no offset * to test, add correction XY * fixed * final fixes * move function out * non-negative scale and dimensions * more optimized renderers, fixed "correction" for rasters * remove specklepy call in helpers * Revert "remove specklepy call in helpers" This reverts commit de002e07f095db88b284cfc0ca4fddad8fa09966. * rely on string type for dataStorage.currentUnits * fix units format * standardize WKT format (#197) * updato to m1 resource on CI, aligned with https://github.com/specklesystems/speckle-sharp/pull/3449 (no rosetta) (#204) * raster renderer edge cases (#206) * simplify and triangulate large meshes (#201) * simplify and triangulate large meshes * more clear logic with adding points * apply coeff to inner rings * fix indexError * max pts * updated links * tags * remove extra readme * syntax error --- .circleci/config.yml | 7 +- metadata.txt | 7 +- plugin_utils/helpers.py | 47 +- speckle/converter/features/GisFeature.py | 13 + .../converter/features/feature_conversions.py | 1598 ++++++++++------- speckle/converter/geometry/conversions.py | 9 +- speckle/converter/geometry/mesh.py | 11 +- speckle/converter/geometry/polygon.py | 19 +- speckle/converter/geometry/utils.py | 82 +- speckle/converter/layers/symbology.py | 16 +- speckle/converter/layers/utils.py | 61 +- speckle/utils/project_vars.py | 7 +- speckle_qgis.py | 3 +- .../converter/layer_tests/test_layer_utils.py | 1 - tests/unit/utils/test_validation.py | 7 +- 15 files changed, 1030 insertions(+), 858 deletions(-) create mode 100644 speckle/converter/features/GisFeature.py diff --git a/.circleci/config.yml b/.circleci/config.yml index 30c7fbe2..13bb8398 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -7,6 +7,7 @@ orbs: # Upload artifacts to s3 aws-s3: circleci/aws-s3@2.0.0 codecov: codecov/codecov@3.3.0 + macos: circleci/macos@2.4.1 jobs: @@ -83,9 +84,10 @@ jobs: VERSION=$(echo "$SEMVER" | sed -e 's/[a-zA-Z]*\///') python patch_version.py $VERSION - run: - name: Remove extra License + name: Remove extra License and Readme command: | rm ./specklepy_qt_ui/LICENSE + rm ./specklepy_qt_ui/qt_ui/README.md - run: name: ZIP plugin command: pb_tool zip @@ -213,7 +215,8 @@ jobs: build-installer-mac: macos: - xcode: 12.5.1 + xcode: 13.4.1 + resource_class: macos.m1.medium.gen1 parameters: runtime: type: string diff --git a/metadata.txt b/metadata.txt index bf2e9715..5264557b 100644 --- a/metadata.txt +++ b/metadata.txt @@ -5,6 +5,7 @@ [general] name=Speckle qgisMinimumVersion=3.28.15 +qgisMaximumVersion=3.34.3 description=Speckle 2.0 Connector for QGIS version=0.0.99 author=Speckle Systems @@ -14,11 +15,11 @@ about= The Speckle QGIS plugin allows you send and receive data from multiple sources to/from several layers in your project, and store their geometry (as well as their contained metadata), in a Speckle server. Don't know what Speckle is? You're not alone! Find out more at https://speckle.systems - You can start exploring by checking out this commit generated in QGIS from open web data sources (OSM, Google Satellite Tiles and Mapzen Terrain Tiles): https://speckle.xyz/streams/5feae56049/commits/b5f4cc5f3b + You can start exploring by checking out this commit generated in QGIS from open web data sources (OSM, Google Satellite Tiles and Mapzen Terrain Tiles): https://app.speckle.systems/projects/5feae56049/models/1a95ec93ec Requirements: - Speckle Manager: You can download the installer here -> https://speckle.guide/user/manager - - An account in a Speckle server. If you don't have one, feel free to use our public server https://speckle.xyz + - An account in a Speckle server. If you don't have one, feel free to use our public server https://app.speckle.systems/ - An account added in the Accounts section of the Speckle Manager - Windows and Mac compatible @@ -48,7 +49,7 @@ hasProcessingProvider=no # changelog= # Tags are comma separated with spaces allowed -tags=python, speckle, interoperability, collaboration, bim, cad +tags=python, speckle, interoperability, collaboration, 3d, bim, cad, online, server, web, cloud homepage=http://speckle.systems category=Web diff --git a/plugin_utils/helpers.py b/plugin_utils/helpers.py index cd526dd2..399175e3 100644 --- a/plugin_utils/helpers.py +++ b/plugin_utils/helpers.py @@ -5,6 +5,9 @@ import inspect from difflib import SequenceMatcher +from specklepy.objects.units import get_units_from_string +from specklepy.objects.units import get_scale_factor_to_meters + SYMBOL = "_x_x_" UNSUPPORTED_PROVIDERS = ["WFS", "wms", "wcs", "vectortile"] @@ -31,57 +34,25 @@ def get_scale_factor(units: str, dataStorage) -> float: return scale_to_meter -def get_scale_factor_to_meter(units: str) -> float: +def get_scale_factor_to_meter(units_src: str) -> float: try: - unit_scale = { - "meters": 1.0, - "centimeters": 0.01, - "millimeters": 0.001, - "inches": 0.0254, - "feet": 0.3048, - "kilometers": 1000.0, - "mm": 0.001, - "cm": 0.01, - "m": 1.0, - "km": 1000.0, - "in": 0.0254, - "ft": 0.3048, - "yd": 0.9144, - "mi": 1609.340, - } - if ( - units is not None - and isinstance(units, str) - and units.lower() in unit_scale.keys() - ): - return unit_scale[units] + units = get_units_from_string(units_src) + return get_scale_factor_to_meters(units) + except: try: from speckle.utils.panel_logging import logToUser logToUser( - f"Units {units} are not supported. Meters will be applied by default.", + f"Units {units_src} are not supported. Meters will be applied by default.", level=1, func=inspect.stack()[0][3], ) return 1.0 except: print( - f"Units {units} are not supported. Meters will be applied by default." - ) - return 1.0 - except Exception as e: - try: - from speckle.utils.panel_logging import logToUser - - logToUser( - f"{e}. Meters will be applied by default.", - level=2, - func=inspect.stack()[0][3], + f"Units {units_src} are not supported. Meters will be applied by default." ) return 1.0 - except: - print(f"{e}. Meters will be applied by default.") - return 1.0 def jsonFromList(jsonObj: dict, levels: list): diff --git a/speckle/converter/features/GisFeature.py b/speckle/converter/features/GisFeature.py new file mode 100644 index 00000000..0b16b1cb --- /dev/null +++ b/speckle/converter/features/GisFeature.py @@ -0,0 +1,13 @@ +from typing import List, Optional + +from specklepy.objects.base import Base + + +class GisFeature( + Base, speckle_type="Objects.GIS.GisFeature", detachable={"displayValue"} +): + """GIS Feature""" + + geometry: Optional[List[Base]] = None + attributes: Base + displayValue: Optional[List[Base]] = None diff --git a/speckle/converter/features/feature_conversions.py b/speckle/converter/features/feature_conversions.py index f092c24a..a713b7cb 100644 --- a/speckle/converter/features/feature_conversions.py +++ b/speckle/converter/features/feature_conversions.py @@ -2,13 +2,14 @@ import inspect import math import os -from typing import List, Union +from typing import Dict, List, Union import numpy as np import hashlib import scipy as sp -from plugin_utils.helpers import findOrCreatePath, get_scale_factor_to_meter +from plugin_utils.helpers import findOrCreatePath +from speckle.converter.features.GisFeature import GisFeature from speckle.converter.geometry import transform from speckle.converter.geometry.conversions import ( convertToNative, @@ -19,10 +20,8 @@ from speckle.converter.geometry.utils import apply_pt_offsets_rotation_on_send from speckle.converter.layers.utils import ( generate_qgis_app_id, - get_raster_stats, getArrayIndicesFromXY, getElevationLayer, - getHeightWithRemainderFromArray, getRasterArrays, getVariantFromValue, getXYofArrayPoint, @@ -79,15 +78,18 @@ def featureToSpeckle( iterations = 0 try: geom = None + new_geom = None if geomType == "None": - geom = GisNonGeometryElement() + geom = GisNonGeometryElement() # redundant, delete in refactor + new_geom = GisFeature() new_report = {"obj_type": geom.speckle_type, "errors": ""} else: # Try to extract geometry skipped_msg = f"'{geomType}' feature skipped due to invalid geometry" try: geom, iterations = convertToSpeckle(f, selectedLayer, dataStorage) + if geom is not None and geom != "None": if not isinstance(geom.geometry, List): logToUser( @@ -97,6 +99,16 @@ def featureToSpeckle( ) return None + # geom is GisPointElement, GisLineElement, GisPolygonElement + new_geom = GisFeature() + new_geom.geometry = [] + for g in geom.geometry: + obj = g + if isinstance(g, GisPolygonGeometry): + new_geom.displayValue = [] + obj = GisPolygonGeometry(boundary=g.boundary, voids=g.voids) + new_geom.geometry.append(obj) + all_errors = "" for g in geom.geometry: if g is None or g == "None": @@ -114,6 +126,7 @@ def featureToSpeckle( func=inspect.stack()[0][3], ) elif iterations is not None and iterations > 0: + new_geom.displayValue.extend(g.displayValue) all_errors += ( "Polygon display mesh is simplified" + ", " ) @@ -122,6 +135,8 @@ def featureToSpeckle( level=1, func=inspect.stack()[0][3], ) + else: + new_geom.displayValue.extend(g.displayValue) if len(geom.geometry) == 0: all_errors = "No geometry converted" @@ -129,7 +144,7 @@ def featureToSpeckle( {"obj_type": geom.speckle_type, "errors": all_errors} ) - else: # geom is None + else: # geom is None, should not happen, but we should pass the object with attributes anyway new_report = {"obj_type": "", "errors": skipped_msg} logToUser(skipped_msg, level=2, func=inspect.stack()[0][3]) geom = GisNonGeometryElement() @@ -162,11 +177,12 @@ def featureToSpeckle( # if geom is not None and geom!="None": geom.attributes = attributes + new_geom.attributes = attributes dataStorage.latestActionFeaturesReport[ len(dataStorage.latestActionFeaturesReport) - 1 ].update(new_report) - return geom + return new_geom except Exception as e: new_report.update({"errors": e}) @@ -174,7 +190,643 @@ def featureToSpeckle( len(dataStorage.latestActionFeaturesReport) - 1 ].update(new_report) logToUser(e, level=2, func=inspect.stack()[0][3]) - return geom + return new_geom + + +def show_progress(current_row: int, rows: int, layer_name: str, plugin: "SpeckleQGIS"): + """Updates UI with raster conversion %.""" + percentage: int = int(current_row / rows) # i = percentage + if percentage % 10 == 0 or percentage == 5: + logToUser( + f"Converting layer '{layer_name}': {percentage}%...", + level=0, + plugin=plugin.dockwidget, + ) + + +def reproject_raster(layer, crs, resolutionX, resolutionY): + if layer.crs() == crs: + return layer + path = ( + os.path.expandvars(r"%LOCALAPPDATA%") + + "\\Temp\\Speckle_QGIS_temp\\" + + datetime.now().strftime("%Y-%m-%d_%H-%M") + ) + findOrCreatePath(path) + + out = path + "\\out.tif" + file = gdal.Warp( + out, + layer.source(), + dstSRS=crs.authid(), + xRes=resolutionX, + yRes=resolutionY, + ) + return QgsRasterLayer(out, "", "gdal") + + +def get_raster_mesh_coords( + reprojected_raster_stats, + rasterResXY: list, + band1_values: list, + dataStorage, +) -> List[float]: + ( + reprojected_top_right, + reprojectedOriginPt, + reprojectedMaxPt, + reprojected_bottom_left, + rasterResXY_reprojected, + rasterDimensions, + ) = reprojected_raster_stats + + xOrigin = reprojectedOriginPt.x() + yOrigin = reprojectedOriginPt.y() + sizeX = rasterDimensions[0] + x_correction = (reprojected_bottom_left.x() - xOrigin) / sizeX + y_correction = (reprojected_top_right.y() - yOrigin) / rasterDimensions[1] + + list_nested = [ + ( + xOrigin + + rasterResXY[0] + * (ind % sizeX) # ind%sizeX = current item index in the row; +1 = next item + + x_correction * (ind % sizeX), + yOrigin + + rasterResXY[1] + * math.floor( + ind / sizeX + ) # math.floor(ind/sizeX) = current row index; +1 = next row + + y_correction * math.floor(ind / sizeX), + 0, + xOrigin + rasterResXY[0] * (ind % sizeX) + x_correction * (ind % sizeX), + yOrigin + + rasterResXY[1] * math.floor(ind / sizeX + 1) + + y_correction * math.floor(ind / sizeX + 1), + 0, + xOrigin + + rasterResXY[0] * (ind % sizeX + 1) + + x_correction * (ind % sizeX + 1), + yOrigin + + rasterResXY[1] * math.floor(ind / sizeX + 1) + + y_correction * math.floor(ind / sizeX + 1), + 0, + xOrigin + + rasterResXY[0] * (ind % sizeX + 1) + + x_correction * (ind % sizeX + 1), + yOrigin + + rasterResXY[1] * math.floor(ind / sizeX) + + y_correction * math.floor(ind / sizeX), + 0, + ) + for ind, _ in enumerate(band1_values) + ] + list_flattened = [item for sublist in list_nested for item in sublist] + + return list_flattened + + +def apply_offset_rotation_to_vertices_send(vertices: List[float], dataStorage): + for index in range(int(len(vertices) / 3)): + x, y = apply_pt_offsets_rotation_on_send( + vertices[index], vertices[index + 1], dataStorage + ) + vertices[index] = x + vertices[index + 1] = y + + +def get_raster_colors( + layer, + rasterBandVals, + rasterBandNoDataVal, + rasterBandMinVal, + rasterBandMaxVal, + rendererType, + plugin, +): + list_colors = [] + if rendererType == "multibandcolor": + + # mock values for R,G,B channels + vals_red = [0 for _ in rasterBandVals[0]] + vals_green = [0 for _ in rasterBandVals[0]] + vals_blue = [0 for _ in rasterBandVals[0]] + + vals_range_red = 1 + vals_range_green = 1 + vals_range_blue = 1 + + val_min_red = 0 + val_min_green = 0 + val_min_blue = 0 + + val_na_red = None + val_na_green = None + val_na_blue = None + + # get band index for each color channel + bandRed = int(layer.renderer().redBand()) + bandGreen = int(layer.renderer().greenBand()) + bandBlue = int(layer.renderer().blueBand()) + + # assign correct values to R,G,B channels, where available + for band_index in range(len(rasterBandVals)): + # if statements are not exclusive, as QGIS allows to assugn 1 band to several color channels + if band_index + 1 == bandRed: + vals_red = rasterBandVals[band_index] + vals_range_red = ( + rasterBandMaxVal[band_index] - rasterBandMinVal[band_index] + ) + val_min_red = rasterBandMinVal[band_index] + val_na_red = rasterBandNoDataVal[band_index] + if band_index + 1 == bandGreen: + vals_green = rasterBandVals[band_index] + vals_range_green = ( + rasterBandMaxVal[band_index] - rasterBandMinVal[band_index] + ) + val_min_green = rasterBandMinVal[band_index] + val_na_green = rasterBandNoDataVal[band_index] + if band_index + 1 == bandBlue: + vals_blue = rasterBandVals[band_index] + vals_range_blue = ( + rasterBandMaxVal[band_index] - rasterBandMinVal[band_index] + ) + val_min_blue = rasterBandMinVal[band_index] + val_na_blue = rasterBandNoDataVal[band_index] + + list_colors = [ + ( + (255 << 24) + | (int(255 * (vals_red[ind] - val_min_red) / vals_range_red) << 16) + | (int(255 * (vals_green[ind] - val_min_green) / vals_range_green) << 8) + | int(255 * (vals_blue[ind] - val_min_blue) / vals_range_blue) + if ( + vals_red[ind] != val_na_red + and vals_green[ind] != val_na_green + and vals_blue[ind] != val_na_blue + ) + else (0 << 24) + (0 << 16) + (0 << 8) + 0 + ) + for ind, _ in enumerate(rasterBandVals[0]) + for _ in range(4) + ] + + elif rendererType == "paletted": + try: + bandIndex = layer.renderer().band() - 1 # int + renderer_classes = layer.renderer().classes() + class_rgbs = [ + renderer_classes[class_ind].color.getRgb() + for class_ind in range(len(renderer_classes)) + ] + + for val in rasterBandVals[bandIndex]: + rgb = None + color = (0 << 24) + (0 << 16) + (0 << 8) + 0 + + for class_ind in range(len(renderer_classes)): + if rasterBandVals[bandIndex] == rasterBandNoDataVal[bandIndex]: + rgb = None + elif class_ind < len(renderer_classes) - 1: + if val >= float( + renderer_classes[class_ind].value + ) and val < float(renderer_classes[class_ind + 1].value): + rgb = class_rgbs[class_ind] + elif class_ind == len(renderer_classes) - 1 and val >= float( + renderer_classes[class_ind].value + ): + rgb = class_rgbs[class_ind] + + if rgb is not None: + color = (255 << 24) + (rgb[0] << 16) + (rgb[1] << 8) + rgb[2] + + # add color value after looping through classes/categories + list_colors.extend([color, color, color, color]) + + except Exception as e: + # log warning, but don't prevent conversion + logToUser( + f"Raster renderer of type '{rendererType}' couldn't read the renderer class values, default renderer will be applied: {e}", + level=1, + func=inspect.stack()[0][3], + plugin=plugin.dockwidget, + ) + + return get_raster_colors( + layer, + rasterBandVals, + rasterBandNoDataVal, + rasterBandMinVal, + rasterBandMaxVal, + None, + plugin, + ) + + elif rendererType == "singlebandpseudocolor": + try: + bandIndex = layer.renderer().band() - 1 # int + + renderer_classes = layer.renderer().legendSymbologyItems() + class_rgbs = [ + renderer_classes[class_ind][1].getRgb() + for class_ind in range(len(renderer_classes)) + ] + + for val in rasterBandVals[bandIndex]: + rgb = None + color = (0 << 24) + (0 << 16) + (0 << 8) + 0 + + for class_ind in range(len(renderer_classes)): + if rasterBandVals[bandIndex] == rasterBandNoDataVal[bandIndex]: + rgb = None + elif class_ind < len(renderer_classes) - 1: + if val >= float(renderer_classes[class_ind][0]) and val < float( + renderer_classes[class_ind + 1][0] + ): + rgb = class_rgbs[class_ind] + elif class_ind == len(renderer_classes) - 1 and val >= float( + renderer_classes[class_ind][0] + ): + rgb = class_rgbs[class_ind] + + if rgb is not None: + color = (255 << 24) + (rgb[0] << 16) + (rgb[1] << 8) + rgb[2] + + # add color value after looping through classes/categories + list_colors.extend([color, color, color, color]) + + except Exception as e: + # log warning, but don't prevent conversion + logToUser( + f"Raster renderer of type '{rendererType}' couldn't read the renderer class values, default renderer will be applied: {e}", + level=1, + func=inspect.stack()[0][3], + plugin=plugin.dockwidget, + ) + + return get_raster_colors( + layer, + rasterBandVals, + rasterBandNoDataVal, + rasterBandMinVal, + rasterBandMaxVal, + None, + plugin, + ) + + else: # greyscale + if rendererType != "singlebandgray": + logToUser( + f"Raster renderer type {rendererType} is not supported, default renderer will be applied", + level=1, + func=inspect.stack()[0][3], + ) + pixMin = rasterBandMinVal[0] + pixMax = rasterBandMaxVal[0] + vals_range = pixMax - pixMin + + list_colors: List[int] = [ + ( + (255 << 24) + | (int(255 * (rasterBandVals[0][ind] - pixMin) / vals_range) << 16) + | (int(255 * (rasterBandVals[0][ind] - pixMin) / vals_range) << 8) + | int(255 * (rasterBandVals[0][ind] - pixMin) / vals_range) + if rasterBandVals[0][ind] != rasterBandNoDataVal[0] + else (0 << 24) + (0 << 16) + (0 << 8) + 0 + ) + for ind, _ in enumerate(rasterBandVals[0]) + for _ in range(4) + ] + + return list_colors + + +def get_vertices_height( + vertices_filtered: List[float], + xy_z_values: dict, + vertices_list_index: int, + height_array, + indices: tuple, +): + index1, index1_0, index2, index2_0 = indices + # top vertices ###################################### + try: + z1 = xy_z_values[ + ( + vertices_filtered[vertices_list_index], + vertices_filtered[vertices_list_index + 1], + ) + ] + except KeyError: + if index1 > 0 and index2 > 0: + z1 = height_array[index1_0][index2_0] + elif index1 > 0: + z1 = height_array[index1_0][index2] + elif index2 > 0: + z1 = height_array[index1][index2_0] + else: + z1 = height_array[index1][index2] + + if z1 is not None: + xy_z_values[ + ( + vertices_filtered[vertices_list_index], + vertices_filtered[vertices_list_index + 1], + ) + ] = z1 + + #################### z4 + try: + z4 = xy_z_values[ + ( + vertices_filtered[vertices_list_index + 9], + vertices_filtered[vertices_list_index + 10], + ) + ] + except: + if index1 > 0: + z4 = height_array[index1_0][index2] + else: + z4 = height_array[index1][index2] + + if z4 is not None: + xy_z_values[ + ( + vertices_filtered[vertices_list_index + 9], + vertices_filtered[vertices_list_index + 10], + ) + ] = z4 + + # bottom vertices ###################################### + z3 = height_array[index1][index2] + if z3 is not None: + xy_z_values[ + ( + vertices_filtered[vertices_list_index + 6], + vertices_filtered[vertices_list_index + 7], + ) + ] = z3 + + try: + z2 = xy_z_values[ + ( + vertices_filtered[vertices_list_index + 3], + vertices_filtered[vertices_list_index + 4], + ) + ] + except: + if index2 > 0: + z2 = height_array[index1][index2_0] + else: + z2 = height_array[index1][index2] + if z2 is not None: + xy_z_values[ + ( + vertices_filtered[vertices_list_index + 3], + vertices_filtered[vertices_list_index + 4], + ) + ] = z2 + return z1, z2, z3, z4 + + +def get_raster_band_data( + selectedLayer, + ds, + index, + rasterBandNames, + rasterBandNoDataVal, + rasterBandVals, + rasterBandMinVal, + rasterBandMaxVal, +) -> List[float]: + rasterBandNames.append(selectedLayer.bandName(index + 1)) + rb = ds.GetRasterBand(index + 1) + valMin = ( + selectedLayer.dataProvider() + .bandStatistics(index + 1, QgsRasterBandStats.All) + .minimumValue + ) + valMax = ( + selectedLayer.dataProvider() + .bandStatistics(index + 1, QgsRasterBandStats.All) + .maximumValue + ) + bandVals = rb.ReadAsArray().tolist() + + bandValsFlat = [] + [bandValsFlat.extend(item) for item in bandVals] + # look at mesh chunking + + const = float(-1 * math.pow(10, 30)) + defaultNoData = rb.GetNoDataValue() + # print(type(rb.GetNoDataValue())) + + # check whether NA value is too small or raster has too small values + # assign min value of an actual list; re-assign NA val; replace list items to new NA val + try: + # create "safe" fake NA value; replace extreme values with it + fakeNA = max(bandValsFlat) + 1 + bandValsFlatFake = [ + fakeNA if val <= const else val for val in bandValsFlat + ] # replace all values corresponding to NoData value + + # if default NA value is too small + if ( + isinstance(defaultNoData, float) or isinstance(defaultNoData, int) + ) and defaultNoData < const: + # find and rewrite min of actual band values; create new NA value + valMin = min(bandValsFlatFake) + noDataValNew = valMin - 1000 # use new adequate value + rasterBandNoDataVal.append(noDataValNew) + # replace fake NA with new NA + bandValsFlat = [ + noDataValNew if val == fakeNA else val for val in bandValsFlatFake + ] # replace all values corresponding to NoData value + + # if default val unaccessible and minimum val is too small + elif ( + isinstance(defaultNoData, str) or defaultNoData is None + ) and valMin < const: # if there are extremely small values but default NA unaccessible + noDataValNew = valMin + rasterBandNoDataVal.append(noDataValNew) + # replace fake NA with new NA + bandValsFlat = [ + noDataValNew if val == fakeNA else val for val in bandValsFlatFake + ] # replace all values corresponding to NoData value + # last, change minValto actual one + valMin = min(bandValsFlatFake) + + else: + rasterBandNoDataVal.append(rb.GetNoDataValue()) + except: + rasterBandNoDataVal.append(rb.GetNoDataValue()) + + rasterBandVals.append(bandValsFlat) + rasterBandMinVal.append(valMin) + rasterBandMaxVal.append(valMax) + return bandValsFlat + + +def get_height_array_from_elevation_layer(elevationLayer): + elevation_arrays, _, _, all_na = getRasterArrays(elevationLayer) + if elevation_arrays is None: + return None + array_band = elevation_arrays[0] + + const = float(-1 * math.pow(10, 30)) + height_array = np.where( + (array_band < const) | (array_band > -1 * const) | (array_band == all_na[0]), + np.nan, + array_band, + ) + try: + height_array = height_array.astype(float) + except: + try: + arr = [] + for row in height_array: + new_row = [] + for item in row: + try: + new_row.append(float(item)) + except: + new_row.append(np.nan) + arr.append(new_row) + height_array = np.array(arr).astype(float) + except: + height_array = height_array[[isinstance(i, float) for i in height_array]] + return height_array + + +def get_raster_reprojected_stats( + project, projectCRS, selectedLayer, originX, originY, rasterResXY, rasterDimensions +): + # get 4 corners of raster + raster_top_right = QgsPointXY( + originX + rasterResXY[0] * rasterDimensions[0], originY + ) + rasterOriginPoint = QgsPointXY(originX, originY) + rasterMaxPt = QgsPointXY( + originX + rasterResXY[0] * rasterDimensions[0], + originY + rasterResXY[1] * rasterDimensions[1], + ) + raster_bottom_left = QgsPointXY( + originX, + originY + rasterResXY[1] * rasterDimensions[1], + ) + # reproject corners to the project CRS + scale_factor_x = scale_factor_y = 1 + if selectedLayer.crs() != projectCRS: + reprojected_top_right = transform.transform( + project, raster_top_right, selectedLayer.crs(), projectCRS + ) + reprojectedOriginPt = transform.transform( + project, rasterOriginPoint, selectedLayer.crs(), projectCRS + ) + reprojectedMaxPt = transform.transform( + project, rasterMaxPt, selectedLayer.crs(), projectCRS + ) + reprojected_bottom_left = transform.transform( + project, raster_bottom_left, selectedLayer.crs(), projectCRS + ) + + scale_factor_x = abs( + (reprojected_top_right.x() - reprojectedOriginPt.x()) + / (raster_top_right.x() - rasterOriginPoint.x()) + ) + scale_factor_y = abs( + (reprojected_top_right.y() - reprojectedMaxPt.y()) + / (raster_top_right.y() - rasterMaxPt.y()) + ) + + rasterResXY_reprojected = [ + rasterResXY[0] * scale_factor_x, + rasterResXY[1] * scale_factor_y, + ] + rasterDimensions_reprojected = ( + abs( + round( + (reprojected_top_right.x() - reprojectedOriginPt.x()) + / rasterResXY_reprojected[0] + ) + ), + abs( + round( + (reprojected_top_right.y() - reprojectedMaxPt.y()) + / rasterResXY_reprojected[1] + ) + ), + ) + return ( + reprojected_top_right, + reprojectedOriginPt, + reprojectedMaxPt, + reprojected_bottom_left, + rasterResXY_reprojected, + rasterDimensions_reprojected, + ) + + +def get_elevation_indices( + texture_transform, + rasterResXY, + rasterResXY_reprojected, + reprojectedOriginX, + reprojectedOriginY, + h, + v, + elevationResX, + elevationResY, + elevationOriginX, + elevationOriginY, + elevationSizeX, + elevationSizeY, +): + if texture_transform is True: # texture + # index1: index on y-scale + posX, posY = getXYofArrayPoint( + rasterResXY_reprojected, + reprojectedOriginX, + reprojectedOriginY, + h, + v, + ) + + index1, index2, _, _ = getArrayIndicesFromXY( + ( + elevationResX, + elevationResY, + elevationOriginX, + elevationOriginY, + elevationSizeX, + elevationSizeY, + None, + None, + ), + posX, + posY, + ) + + index1_0, index2_0, _, _ = getArrayIndicesFromXY( + ( + elevationResX, + elevationResY, + elevationOriginX, + elevationOriginY, + elevationSizeX, + elevationSizeY, + None, + None, + ), + posX - rasterResXY[0], + posY - rasterResXY[1], + ) + else: # elevation + index1 = v + index1_0 = v - 1 + index2 = h + index2_0 = h - 1 + return index1, index1_0, index2, index2_0 def rasterFeatureToSpeckle( @@ -189,9 +841,8 @@ def rasterFeatureToSpeckle( b = GisRasterElement(units=dataStorage.currentUnits) try: - terrain_transform = False - texture_transform = False - # height_list = rasterBandVals[0] + time0 = datetime.now() + terrain_transform = isAppliedLayerTransformByKeywords( selectedLayer, ["elevation", "mesh"], ["texture"], dataStorage ) @@ -202,111 +853,53 @@ def rasterFeatureToSpeckle( b = GisTopography(units=dataStorage.currentUnits) rasterBandCount = selectedLayer.bandCount() - rasterBandNames = [] rasterDimensions = [selectedLayer.width(), selectedLayer.height()] - # if rasterDimensions[0]*rasterDimensions[1] > 1000000 : - # logToUser("Large layer: ", level = 1, func = inspect.stack()[0][3]) ds = gdal.Open(selectedLayer.source(), gdal.GA_ReadOnly) - if ds is None: - return None + rasterResXY = [float(ds.GetGeoTransform()[1]), float(ds.GetGeoTransform()[5])] originX = ds.GetGeoTransform()[0] originY = ds.GetGeoTransform()[3] - rasterOriginPoint = QgsPointXY(originX, originY) - rasterResXY = [float(ds.GetGeoTransform()[1]), float(ds.GetGeoTransform()[5])] - rasterWkt = ds.GetProjection() - rasterProj = ( - QgsCoordinateReferenceSystem.fromWkt(rasterWkt) - .toProj() - .replace(" +type=crs", "") + + reprojected_raster_stats = get_raster_reprojected_stats( + project, + projectCRS, + selectedLayer, + originX, + originY, + rasterResXY, + rasterDimensions, ) + ( + reprojected_top_right, + reprojectedOriginPt, + reprojectedMaxPt, + reprojected_bottom_left, + rasterResXY_reprojected, + rasterDimensions_reprojected, + ) = reprojected_raster_stats + reprojectedOriginX, reprojectedOriginY = ( + reprojectedOriginPt.x(), + reprojectedOriginPt.y(), + ) + + # fill band values rasterBandNoDataVal = [] rasterBandMinVal = [] rasterBandMaxVal = [] rasterBandVals = [] - - # Try to extract geometry - reprojectedPt = QgsGeometry.fromPointXY(QgsPointXY()) - try: - reprojectedPt = rasterOriginPoint - if selectedLayer.crs() != projectCRS: - reprojectedPt = transform.transform( - project, rasterOriginPoint, selectedLayer.crs(), projectCRS - ) - except Exception as error: - # logToUser("Error converting point geometry: " + str(error), level = 2, func = inspect.stack()[0][3]) - logToUser("Error converting point geometry: " + str(error), level=2) - + rasterBandNames = [] for index in range(rasterBandCount): - rasterBandNames.append(selectedLayer.bandName(index + 1)) - rb = ds.GetRasterBand(index + 1) - valMin = ( - selectedLayer.dataProvider() - .bandStatistics(index + 1, QgsRasterBandStats.All) - .minimumValue + bandValsFlat = get_raster_band_data( + selectedLayer, + ds, + index, + rasterBandNames, + rasterBandNoDataVal, + rasterBandVals, + rasterBandMinVal, + rasterBandMaxVal, ) - valMax = ( - selectedLayer.dataProvider() - .bandStatistics(index + 1, QgsRasterBandStats.All) - .maximumValue - ) - bandVals = rb.ReadAsArray().tolist() - - bandValsFlat = [] - [bandValsFlat.extend(item) for item in bandVals] - # look at mesh chunking - - const = float(-1 * math.pow(10, 30)) - defaultNoData = rb.GetNoDataValue() - # print(type(rb.GetNoDataValue())) - - # check whether NA value is too small or raster has too small values - # assign min value of an actual list; re-assign NA val; replace list items to new NA val - try: - # create "safe" fake NA value; replace extreme values with it - fakeNA = max(bandValsFlat) + 1 - bandValsFlatFake = [ - fakeNA if val <= const else val for val in bandValsFlat - ] # replace all values corresponding to NoData value - - # if default NA value is too small - if ( - isinstance(defaultNoData, float) or isinstance(defaultNoData, int) - ) and defaultNoData < const: - # find and rewrite min of actual band values; create new NA value - valMin = min(bandValsFlatFake) - noDataValNew = valMin - 1000 # use new adequate value - rasterBandNoDataVal.append(noDataValNew) - # replace fake NA with new NA - bandValsFlat = [ - noDataValNew if val == fakeNA else val - for val in bandValsFlatFake - ] # replace all values corresponding to NoData value - - # if default val unaccessible and minimum val is too small - elif ( - isinstance(defaultNoData, str) or defaultNoData is None - ) and valMin < const: # if there are extremely small values but default NA unaccessible - noDataValNew = valMin - rasterBandNoDataVal.append(noDataValNew) - # replace fake NA with new NA - bandValsFlat = [ - noDataValNew if val == fakeNA else val - for val in bandValsFlatFake - ] # replace all values corresponding to NoData value - # last, change minValto actual one - valMin = min(bandValsFlatFake) - - else: - rasterBandNoDataVal.append(rb.GetNoDataValue()) - - except: - rasterBandNoDataVal.append(rb.GetNoDataValue()) - - rasterBandVals.append(bandValsFlat) - rasterBandMinVal.append(valMin) - rasterBandMaxVal.append(valMax) b["@(10000)" + selectedLayer.bandName(index + 1) + "_values"] = ( bandValsFlat # [0:int(max_values/rasterBandCount)] ) @@ -316,139 +909,108 @@ def rasterFeatureToSpeckle( b.x_size = rasterDimensions[0] b.y_size = rasterDimensions[1] b.x_origin, b.y_origin = apply_pt_offsets_rotation_on_send( - reprojectedPt.x(), reprojectedPt.y(), dataStorage + reprojectedOriginPt.x(), reprojectedOriginPt.y(), dataStorage ) b.band_count = rasterBandCount b.band_names = rasterBandNames b.noDataValue = rasterBandNoDataVal - # creating a mesh - count = 0 - rendererType = selectedLayer.renderer().type() - - xy_list = [] - z_list = [] - # print(rendererType) - # identify symbology type and if Multiband, which band is which color + # creating a mesh + xy_z_values: Dict[tuple, float] = {} ############################################################# elevationLayer = None - elevationProj = None - if texture_transform is True: - elevationLayer = getElevationLayer(dataStorage) - elif terrain_transform is True: + height_array = None + + if terrain_transform is True: # same layer, copy props elevationLayer = selectedLayer + elevationResX, elevationResY = rasterResXY_reprojected + elevationOriginX, elevationOriginY = ( + reprojectedOriginPt.x(), + reprojectedOriginPt.y(), + ) + elevationSizeX, elevationSizeY = ( + rasterDimensions_reprojected[0], + rasterDimensions_reprojected[1], + ) - if elevationLayer is not None: - settings_elevation_layer = get_raster_stats(elevationLayer) - ( - elevationResX, - elevationResY, - elevationOriginX, - elevationOriginY, - elevationSizeX, - elevationSizeY, - elevationWkt, - elevationProj, - ) = settings_elevation_layer + elif texture_transform is True: + elevation_layer_original: "QgsRasterLayer" = getElevationLayer(dataStorage) - # reproject the elevation layer - if ( - elevationProj is not None - and rasterProj is not None - and elevationProj != rasterProj - ): - try: - p = ( - os.path.expandvars(r"%LOCALAPPDATA%") - + "\\Temp\\Speckle_QGIS_temp\\" - + datetime.now().strftime("%Y-%m-%d_%H-%M") - ) - findOrCreatePath(p) - path = p - out = p + "\\out.tif" - gdal.Warp( - out, - elevationLayer.source(), - dstSRS=selectedLayer.crs().authid(), - xRes=elevationResX, - yRes=elevationResY, - ) + if elevation_layer_original is None: + elevationResX = elevationResY = elevationOriginX = elevationOriginY = ( + elevationSizeX + ) = elevationSizeY = None - elevationLayer = QgsRasterLayer(out, "", "gdal") - settings_elevation_layer = get_raster_stats(elevationLayer) - ( - elevationResX, - elevationResY, - elevationOriginX, - elevationOriginY, - elevationSizeX, - elevationSizeY, - elevationWkt, - elevationProj, - ) = settings_elevation_layer - except Exception as e: - logToUser(f"Reprojection did not succeed: {e}", level=0) - elevation_arrays, all_mins, all_maxs, all_na = getRasterArrays( - elevationLayer - ) - array_band = elevation_arrays[0] - - height_array = np.where( - (array_band < const) - | (array_band > -1 * const) - | (array_band == all_na[0]), - np.nan, - array_band, - ) - try: - height_array = height_array.astype(float) - except: - try: - arr = [] - for row in height_array: - new_row = [] - for item in row: - try: - new_row.append(float(item)) - except: - new_row.append(np.nan) - arr.append(new_row) - height_array = np.array(arr).astype(float) - except: - height_array = height_array[ - [isinstance(i, float) for i in height_array] - ] - else: - elevation_arrays = all_mins = all_maxs = all_na = None - elevationResX = elevationResY = elevationOriginX = elevationOriginY = ( - elevationSizeX - ) = elevationSizeY = elevationWkt = None - height_array = None + logToUser( + f"Elevation layer is not found. Texture transformation for layer '{selectedLayer.name()}' will not be applied", + level=1, + plugin=plugin.dockwidget, + ) + else: + # match elevation layer props to the reprojected SelectedLayer + ds_elevation_original = gdal.Open( + elevation_layer_original.source(), gdal.GA_ReadOnly + ) + + elevation_original_dimensions = [ + elevation_layer_original.width(), + elevation_layer_original.height(), + ] + elevation_original_ResXY = [ + float(ds_elevation_original.GetGeoTransform()[1]), + float(ds_elevation_original.GetGeoTransform()[5]), + ] + elevation_original_originX = ds_elevation_original.GetGeoTransform()[0] + elevation_original_originY = ds_elevation_original.GetGeoTransform()[3] + + reprojected_elevation_stats = get_raster_reprojected_stats( + project, + projectCRS, + elevation_layer_original, + elevation_original_originX, + elevation_original_originY, + elevation_original_ResXY, + elevation_original_dimensions, + ) + + ( + elevation_reprojected_top_right, + elevation_reprojectedOriginPt, + elevation_reprojectedMaxPt, + elevation_reprojected_bottom_left, + elevation_resXY_reprojected, + elevation_dimensions_reprojected, + ) = reprojected_elevation_stats + + elevationOriginX = elevation_reprojectedOriginPt.x() + elevationOriginY = elevation_reprojectedOriginPt.y() + + # overwrite resolution & dimension to match raster + elevationResX, elevationResY = ( + elevation_resXY_reprojected[0], + elevation_resXY_reprojected[1], + ) + elevationSizeX, elevationSizeY = ( + elevation_dimensions_reprojected[0], + elevation_dimensions_reprojected[1], + ) + elevationLayer = elevation_layer_original + ################# + + if elevationLayer is not None: + height_array = get_height_array_from_elevation_layer(elevationLayer) + if height_array is None: + logToUser( + f"Elevation layer is not found. Texture transformation for layer '{selectedLayer.name()}' will not be applied", + level=1, + plugin=plugin.dockwidget, + ) largeTransform = False - if texture_transform is True and elevationLayer is None: - logToUser( - f"Elevation layer is not found. Texture transformation for layer '{selectedLayer.name()}' will not be applied", - level=1, - plugin=plugin.dockwidget, - ) - elif ( - texture_transform is True - and rasterDimensions[1] * rasterDimensions[0] >= 10000 - and elevationProj is not None - and rasterProj is not None - and elevationProj != rasterProj - ): - # warning if >= 100x100 raster is being projected to an elevation with different CRS - logToUser( - f"Texture transformation for the layer '{selectedLayer.name()}' might take a while 🕒\nTip: reproject one of the layers (texture or elevation) to the other layer's CRS. When both layers have the same CRS, texture transformation will be much faster.", - level=0, - plugin=plugin.dockwidget, - ) - largeTransform = True - elif ( + if ( texture_transform is True + and height_array is not None and rasterDimensions[1] * rasterDimensions[0] >= 250000 ): # warning if >= 500x500 raster is being projected to any elevation @@ -459,440 +1021,108 @@ def rasterFeatureToSpeckle( ) largeTransform = True ############################################################ - faces_array = [] - colors_array = [] - vertices_array = [] - array_z = [] # size is large by 1 than the raster size, in both dimensions - time0 = datetime.now() - for v in range(rasterDimensions[1]): # each row, Y - if largeTransform is True: - if v == int(rasterDimensions[1] / 20): - logToUser( - f"Converting layer '{selectedLayer.name()}': 5%...", - level=0, - plugin=plugin.dockwidget, - ) - elif v == int(rasterDimensions[1] / 10): - logToUser( - f"Converting layer '{selectedLayer.name()}': 10%...", - level=0, - plugin=plugin.dockwidget, - ) - elif v == int(rasterDimensions[1] / 5): - logToUser( - f"Converting layer '{selectedLayer.name()}': 20%...", - level=0, - plugin=plugin.dockwidget, - ) - elif v == int(rasterDimensions[1] * 2 / 5): - logToUser( - f"Converting layer '{selectedLayer.name()}': 40%...", - level=0, - plugin=plugin.dockwidget, - ) - elif v == int(rasterDimensions[1] * 3 / 5): - logToUser( - f"Converting layer '{selectedLayer.name()}': 60%...", - level=0, - plugin=plugin.dockwidget, - ) - elif v == int(rasterDimensions[1] * 4 / 5): - logToUser( - f"Converting layer '{selectedLayer.name()}': 80%...", - level=0, - plugin=plugin.dockwidget, - ) - elif v == int(rasterDimensions[1] * 9 / 10): - logToUser( - f"Converting layer '{selectedLayer.name()}': 90%...", - level=0, - plugin=plugin.dockwidget, - ) - vertices = [] - faces = [] - colors = [] - row_z = [] - row_z_bottom = [] - for h in range(rasterDimensions[0]): # item in a row, X - pt1 = QgsPointXY( - rasterOriginPoint.x() + h * rasterResXY[0], - rasterOriginPoint.y() + v * rasterResXY[1], - ) - pt2 = QgsPointXY( - rasterOriginPoint.x() + h * rasterResXY[0], - rasterOriginPoint.y() + (v + 1) * rasterResXY[1], - ) - pt3 = QgsPointXY( - rasterOriginPoint.x() + (h + 1) * rasterResXY[0], - rasterOriginPoint.y() + (v + 1) * rasterResXY[1], - ) - pt4 = QgsPointXY( - rasterOriginPoint.x() + (h + 1) * rasterResXY[0], - rasterOriginPoint.y() + v * rasterResXY[1], - ) - # first, get point coordinates with correct position and resolution, then reproject each: - if selectedLayer.crs() != projectCRS: - pt1 = transform.transform( - project, src=pt1, crsSrc=selectedLayer.crs(), crsDest=projectCRS - ) - pt2 = transform.transform( - project, src=pt2, crsSrc=selectedLayer.crs(), crsDest=projectCRS - ) - pt3 = transform.transform( - project, src=pt3, crsSrc=selectedLayer.crs(), crsDest=projectCRS - ) - pt4 = transform.transform( - project, src=pt4, crsSrc=selectedLayer.crs(), crsDest=projectCRS - ) + array_z = [] # size is larger by 1 than the raster size, in both dimensions + + faces_filtered = [] + colors_filtered = [] + vertices_filtered = [] - z1 = z2 = z3 = z4 = 0 - index1 = index1_0 = None + # construct mesh + band1_values = rasterBandVals[0] + list_nested = [ + (4, 4 * ind, 4 * ind + 1, 4 * ind + 2, 4 * ind + 3) + for ind, _ in enumerate(band1_values) + ] + faces_filtered = [item for sublist in list_nested for item in sublist] + vertices_filtered = get_raster_mesh_coords( + reprojected_raster_stats, + rasterResXY_reprojected, + band1_values, + dataStorage, + ) + rendererType = selectedLayer.renderer().type() + colors_filtered = get_raster_colors( + selectedLayer, + rasterBandVals, + rasterBandNoDataVal, + rasterBandMinVal, + rasterBandMaxVal, + rendererType, + plugin, + ) + ############################################################################### + if texture_transform is True or terrain_transform is True: + for v in range(rasterDimensions[1]): # each row, Y + if largeTransform is True: + show_progress(v, rasterDimensions[1], selectedLayer.name(), plugin) - ############################################################# - if ( - terrain_transform is True or texture_transform is True - ) and height_array is not None: - if texture_transform is True: # texture - # index1: index on y-scale - posX, posY = getXYofArrayPoint( - ( - rasterResXY[0], - rasterResXY[1], - originX, - originY, - ), + row_z = [] + row_z_bottom = [] + for h in range(rasterDimensions[0]): # item in a row, X + vertices_list_index = 3 * 4 * (v * rasterDimensions[0] + h) + colors_list_index = 4 * (v * rasterDimensions[0] + h) + + z1 = z2 = z3 = z4 = 0 + index1 = index1_0 = None + + ############################################################# + if height_array is not None: + index1, index1_0, index2, index2_0 = get_elevation_indices( + texture_transform, + rasterResXY, + rasterResXY_reprojected, + reprojectedOriginX, + reprojectedOriginY, h, v, - selectedLayer, - elevationLayer, - dataStorage, - ) - - index1, index2, remainder1, remainder2 = getArrayIndicesFromXY( - ( - elevationResX, - elevationResY, - elevationOriginX, - elevationOriginY, - elevationSizeX, - elevationSizeY, - elevationWkt, - elevationProj, - ), - posX, - posY, - ) - ( - index1_0, - index2_0, - remainder1_0, - remainder2_0, - ) = getArrayIndicesFromXY( - ( - elevationResX, - elevationResY, - elevationOriginX, - elevationOriginY, - elevationSizeX, - elevationSizeY, - elevationWkt, - elevationProj, - ), - posX - rasterResXY[0], - posY - rasterResXY[1], + elevationResX, + elevationResY, + elevationOriginX, + elevationOriginY, + elevationSizeX, + elevationSizeY, ) - else: # elevation - index1 = v - index1_0 = v - 1 - index2 = h - index2_0 = h - 1 - - if index1 is None or index1_0 is None: - # count += 4 - # continue # skip the pixel - z1 = z2 = z3 = z4 = np.nan - else: - # top vertices ###################################### - try: - z1 = z_list[xy_list.index((pt1.x(), pt1.y()))] - except: - if index1 > 0 and index2 > 0: - z1 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1_0, index2_0 - ) - elif index1 > 0: - z1 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1_0, index2 - ) - elif index2 > 0: - z1 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1, index2_0 - ) - else: - z1 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1, index2 - ) - - if z1 is not None: - z_list.append(z1) - xy_list.append((pt1.x(), pt1.y())) - - #################### z4 - try: - z4 = z_list[xy_list.index((pt4.x(), pt4.y()))] - except: - if index1 > 0: - z4 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1_0, index2 - ) - else: - z4 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1, index2 - ) - - if z4 is not None: - z_list.append(z4) - xy_list.append((pt4.x(), pt4.y())) - - # bottom vertices ###################################### - z3 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1, index2 - ) - if z3 is not None: - z_list.append(z3) - xy_list.append((pt3.x(), pt3.y())) - - try: - z2 = z_list[xy_list.index((pt2.x(), pt2.y()))] - except: - if index2 > 0: - z2 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1, index2_0 - ) - else: - z2 = getHeightWithRemainderFromArray( - height_array, texture_transform, index1, index2 - ) - if z2 is not None: - z_list.append(z2) - xy_list.append((pt2.x(), pt2.y())) - - ############################################## - - max_len = rasterDimensions[0] * 4 + 4 - if len(z_list) > max_len: - z_list = z_list[len(z_list) - max_len :] - xy_list = xy_list[len(xy_list) - max_len :] - - ### list to smoothen later: - if h == 0: - row_z.append(z1) - row_z_bottom.append(z2) - row_z.append(z4) - row_z_bottom.append(z3) - - ######################################################## - x1, y1 = apply_pt_offsets_rotation_on_send( - pt1.x(), pt1.y(), dataStorage - ) - x2, y2 = apply_pt_offsets_rotation_on_send( - pt2.x(), pt2.y(), dataStorage - ) - x3, y3 = apply_pt_offsets_rotation_on_send( - pt3.x(), pt3.y(), dataStorage - ) - x4, y4 = apply_pt_offsets_rotation_on_send( - pt4.x(), pt4.y(), dataStorage - ) - - vertices.append( - [x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4] - ) ## add 4 points - current_vertices = ( - v * rasterDimensions[0] * 4 + h * 4 - ) # len(np.array(faces_array).flatten()) * 4 / 5 - faces.append( - [ - 4, - current_vertices, - current_vertices + 1, - current_vertices + 2, - current_vertices + 3, - ] - ) - - # color vertices according to QGIS renderer - alpha = 255 - color = (alpha << 24) + (0 << 16) + (0 << 8) + 0 - noValColor: tuple[int] = ( - selectedLayer.renderer().nodataColor().getRgb() - ) # RGB or RGBA, 3 or 4 values - - colorLayer = selectedLayer - currentRasterBandCount = rasterBandCount - - if ( - (terrain_transform is True or texture_transform is True) - and height_array is not None - and (index1 is None or index1_0 is None) - ): # transparent color - alpha = 0 - color = (alpha << 24) + (0 << 16) + (0 << 8) + 0 - elif rendererType == "multibandcolor": - valR = 0 - valG = 0 - valB = 0 - bandRed = int(colorLayer.renderer().redBand()) - bandGreen = int(colorLayer.renderer().greenBand()) - bandBlue = int(colorLayer.renderer().blueBand()) - - for k in range(currentRasterBandCount): - valRange = rasterBandMaxVal[k] - rasterBandMinVal[k] - if valRange == 0: - colorVal = 0 - elif ( - rasterBandVals[k][int(count / 4)] == rasterBandNoDataVal[k] - ): - colorVal = 0 - alpha = 0 - # break + if index1 is None or index1_0 is None: + z1 = z2 = z3 = z4 = np.nan else: - colorVal = int( - ( - rasterBandVals[k][int(count / 4)] - - rasterBandMinVal[k] - ) - / valRange - * 255 - ) - - if k + 1 == bandRed: - valR = colorVal - if k + 1 == bandGreen: - valG = colorVal - if k + 1 == bandBlue: - valB = colorVal - - color = (alpha << 24) + (valR << 16) + (valG << 8) + valB - - elif rendererType == "paletted": - bandIndex = colorLayer.renderer().band() - 1 # int - # if textureLayer is not None: - # value = texture_arrays[bandIndex][index1][index2] - # else: - value = rasterBandVals[bandIndex][ - int(count / 4) - ] # find in the list and match with color - - rendererClasses = colorLayer.renderer().classes() - for c in range(len(rendererClasses) - 1): - if ( - value >= rendererClasses[c].value - and value <= rendererClasses[c + 1].value - ): - rgb = rendererClasses[c].color.getRgb() - color = ( - (255 << 24) + (rgb[0] << 16) + (rgb[1] << 8) + rgb[2] - ) - break - if value == rasterBandNoDataVal[bandIndex]: - alpha = 0 - color = (alpha << 24) + (0 << 16) + (0 << 8) + 0 - - elif rendererType == "singlebandpseudocolor": - bandIndex = colorLayer.renderer().band() - 1 # int - # if textureLayer is not None: - # value = texture_arrays[bandIndex][index1][index2] - # else: - value = rasterBandVals[bandIndex][ - int(count / 4) - ] # find in the list and match with color - - rendererClasses = colorLayer.renderer().legendSymbologyItems() - for c in range(len(rendererClasses) - 1): - if value >= float(rendererClasses[c][0]) and value <= float( - rendererClasses[c + 1][0] - ): - rgb = rendererClasses[c][1].getRgb() - color = ( - (255 << 24) + (rgb[0] << 16) + (rgb[1] << 8) + rgb[2] + z1, z2, z3, z4 = get_vertices_height( + vertices_filtered, + xy_z_values, + vertices_list_index, + height_array, + (index1, index1_0, index2, index2_0), ) - break - if value == rasterBandNoDataVal[bandIndex]: - alpha = 0 - color = (alpha << 24) + (0 << 16) + (0 << 8) + 0 - else: - if rendererType == "singlebandgray": - bandIndex = colorLayer.renderer().grayBand() - 1 - elif rendererType == "hillshade": - bandIndex = colorLayer.renderer().band() - 1 - elif rendererType == "contour": - try: - bandIndex = colorLayer.renderer().inputBand() - 1 - except: - try: - bandIndex = colorLayer.renderer().band() - 1 - except: - bandIndex = 0 - else: # e.g. single band data - bandIndex = 0 - - value = rasterBandVals[bandIndex][int(count / 4)] - if ( - rasterBandMinVal[bandIndex] - <= value - <= rasterBandMaxVal[bandIndex] - ): - # REMAP band values to (0,255) range - valRange = ( - rasterBandMaxVal[bandIndex] - rasterBandMinVal[bandIndex] - ) - if valRange == 0: - colorVal = 0 - else: - colorVal = int( - ( - rasterBandVals[bandIndex][int(count / 4)] - - rasterBandMinVal[bandIndex] - ) - / valRange - * 255 - ) - color = ( - (alpha << 24) - + (colorVal << 16) - + (colorVal << 8) - + colorVal + ### height list to smoothen later: + if h == 0: + row_z.append(z1) + row_z_bottom.append(z2) + row_z.append(z4) + row_z_bottom.append(z3) + + # color vertices according to QGIS renderer + if height_array is not None and ( + index1 is None or index1_0 is None + ): # transparent color + colors_filtered[colors_list_index] = colors_filtered[ + colors_list_index + 1 + ] = colors_filtered[colors_list_index + 2] = colors_filtered[ + colors_list_index + 3 + ] = ( + (0 << 24) + (0 << 16) + (0 << 8) + 0 ) - elif value == rasterBandNoDataVal[bandIndex]: - alpha = 0 - color = (alpha << 24) + (0 << 16) + (0 << 8) + 0 - colors.append([color, color, color, color]) - count += 4 + if terrain_transform is True or texture_transform is True: + if v == 0: + array_z.append(row_z) + array_z.append(row_z_bottom) - # after each row - vertices_array.append(vertices) - faces_array.append(faces) - colors_array.append(colors) - - if v == 0: - array_z.append(row_z) - array_z.append(row_z_bottom) - - time1 = datetime.now() - # print(f"Time to get Raster: {(time1-time0).total_seconds()} sec") - # after the entire loop - faces_filtered = [] - colors_filtered = [] - vertices_filtered = [] - - ## end of the the table - smooth = False - if terrain_transform is True or texture_transform is True: - smooth = True - if smooth is True and len(row_z) > 2 and len(array_z) > 2: + ## smoothen z-values + if ( + (terrain_transform is True or texture_transform is True) + and len(row_z) > 2 + and len(array_z) > 2 + ): array_z_nans = np.array(array_z) array_z_filled = np.array(array_z) @@ -905,62 +1135,31 @@ def rasterFeatureToSpeckle( if texture_transform is True: sigma = 1 # for texture - # increase sigma if needed - try: - unitsRaster = QgsUnitTypes.encodeUnit( - selectedLayer.crs().mapUnits() - ) - unitsElevation = QgsUnitTypes.encodeUnit( - elevationLayer.crs().mapUnits() - ) - # print(unitsRaster) - # print(unitsElevation) - resRasterX = get_scale_factor_to_meter(unitsRaster) * rasterResXY[0] - resElevX = get_scale_factor_to_meter(unitsElevation) * elevationResX - # print(resRasterX) - # print(resElevX) - if resRasterX / resElevX >= 2 or resElevX / resRasterX >= 2: - sigma = math.sqrt( - max(resRasterX / resElevX, resElevX / resRasterX) - ) - # print(sigma) - except: - pass - gaussian_array = sp.ndimage.filters.gaussian_filter( array_z_filled, sigma, mode="nearest" ) + # update vertices_filtered with z-value for v in range(rasterDimensions[1]): # each row, Y for h in range(rasterDimensions[0]): # item in a row, X + vertices_list_index = 3 * 4 * (v * rasterDimensions[0] + h) if not np.isnan(array_z_nans[v][h]): - vertices_item = vertices_array[v][h] - # print(vertices_item) - vertices_item[2] = gaussian_array[v][h] - vertices_item[5] = gaussian_array[v + 1][h] - vertices_item[8] = gaussian_array[v + 1][h + 1] - vertices_item[11] = gaussian_array[v][h + 1] - vertices_filtered.extend(vertices_item) - - currentFaces = len(faces_filtered) / 5 * 4 - faces_filtered.extend( - [ - 4, - currentFaces, - currentFaces + 1, - currentFaces + 2, - currentFaces + 3, - ] - ) - # print(faces_filtered) - colors_filtered.extend(colors_array[v][h]) - # print(colors_array[v][h]) - else: - faces_filtered = np.array(faces_array).flatten().tolist() - colors_filtered = np.array(colors_array).flatten().tolist() - vertices_filtered = np.array(vertices_array).flatten().tolist() + vertices_filtered[vertices_list_index + 2] = gaussian_array[v][ + h + ] + vertices_filtered[vertices_list_index + 5] = gaussian_array[ + v + 1 + ][h] + vertices_filtered[vertices_list_index + 8] = gaussian_array[ + v + 1 + ][h + 1] + vertices_filtered[vertices_list_index + 11] = gaussian_array[v][ + h + 1 + ] + + # apply offset & rotation + apply_offset_rotation_to_vertices_send(vertices_filtered, dataStorage) - # if len(colors)/4*5 == len(faces) and len(colors)*3 == len(vertices): mesh = constructMeshFromRaster( vertices_filtered, faces_filtered, colors_filtered, dataStorage ) @@ -974,6 +1173,8 @@ def rasterFeatureToSpeckle( plugin=plugin.dockwidget, ) + time1 = datetime.now() + print(f"Time to convert Raster: {(time1-time0).total_seconds()} sec") return b except Exception as e: @@ -987,7 +1188,9 @@ def featureToNative(feature: Base, fields: "QgsFields", dataStorage): try: qgsGeom = None - if isinstance(feature, GisNonGeometryElement): + if isinstance(feature, GisNonGeometryElement) or ( + isinstance(feature, GisFeature) and feature.geometry is None + ): pass else: try: @@ -1006,7 +1209,14 @@ def featureToNative(feature: Base, fields: "QgsFields", dataStorage): qgsGeom = convertToNative(speckle_geom, dataStorage) elif isinstance(speckle_geom, list): - if len(speckle_geom) == 1: + # add condition for new GisFeature class + if ( + isinstance(feature, GisFeature) + and isinstance(speckle_geom[0], GisPolygonGeometry) + and speckle_geom[0].boundary is None + ): + qgsGeom = convertToNativeMulti(feature.displayValue, dataStorage) + elif len(speckle_geom) == 1: qgsGeom = convertToNative(speckle_geom[0], dataStorage) elif len(speckle_geom) > 1: qgsGeom = convertToNativeMulti(speckle_geom, dataStorage) diff --git a/speckle/converter/geometry/conversions.py b/speckle/converter/geometry/conversions.py index 82a4d761..b8123ae3 100644 --- a/speckle/converter/geometry/conversions.py +++ b/speckle/converter/geometry/conversions.py @@ -237,9 +237,7 @@ def convertToSpeckle( v[1] for v in enumerate(poly.exteriorRing().vertices()) ] except: - boundaryPts = [ - v[1] for v in enumerate(poly.vertices()) - ] + boundaryPts = [v[1] for v in enumerate(poly.vertices())] if height is not None: if isFlat(boundaryPts) is False: logToUser( @@ -445,6 +443,9 @@ def convertToNativeMulti( return multiPolylineToNative(items, dataStorage) # elif isinstance(first, Arc) or isinstance(first, Polycurve) or isinstance(first, Ellipse) or isinstance(first, Circle) or isinstance(first, Curve): # return [convertToNative(it, dataStorage) for it in items] + elif isinstance(first, Mesh): + converted: QgsMultiPolygon = meshToNative(items, dataStorage) + return converted elif isinstance(first, Base): try: displayVals = [] @@ -453,7 +454,7 @@ def convertToNativeMulti( displayVals.extend(it.displayValue) except: displayVals.extend(it["@displayValue"]) - if isinstance(first, GisPolygonGeometry): + if isinstance(first, GisPolygonGeometry) or isinstance(first, Mesh): if first.boundary is None: converted: QgsMultiPolygon = meshToNative( displayVals, dataStorage diff --git a/speckle/converter/geometry/mesh.py b/speckle/converter/geometry/mesh.py index cbde658e..8aafe936 100644 --- a/speckle/converter/geometry/mesh.py +++ b/speckle/converter/geometry/mesh.py @@ -1,4 +1,5 @@ import inspect +import math from typing import List, Tuple, Union from specklepy.objects.geometry import Mesh, Point from specklepy.objects.other import RenderMaterial @@ -212,16 +213,14 @@ def meshPartsFromPolygon( coef = 1 maxPoints = 5000 if len(polyBorder) >= maxPoints: - coef = int(len(polyBorder) / maxPoints) + coef = math.floor(len(polyBorder) / maxPoints) if len(polyBorder) < 3: return None, None, None, None, None col = featureColorfromNativeRenderer(feature, layer) - if ( - len(voidsAsPts) == 0 - ): # only if there is a mesh with no voids and large amount of points + if len(voidsAsPts) == 0: # only if there is a mesh with no voids # floor: need positive - clockwise (looking down); cap need negative (counter-clockwise) polyBorder = fix_orientation(polyBorder, True, coef) # clockwise @@ -229,8 +228,8 @@ def meshPartsFromPolygon( polyBorder.reverse() # when no extrusions: face up, counter-clockwise for k, ptt in enumerate(polyBorder): # pointList: - pt = polyBorder[k * coef] if k < maxPoints: + pt = polyBorder[k * coef] vertices.extend([pt.x, pt.y, pt.z]) total_vertices += 1 else: @@ -351,7 +350,7 @@ def meshPartsFromPolygon( # get points from original geometry ################################# triangulated_geom, vertices3d_original, iterations = triangulatePolygon( - feature_geom, dataStorage, xform + feature_geom, dataStorage, coef, xform ) # temporary solution, as the list of points is not the same anymore: diff --git a/speckle/converter/geometry/polygon.py b/speckle/converter/geometry/polygon.py index 130fb019..c721e489 100644 --- a/speckle/converter/geometry/polygon.py +++ b/speckle/converter/geometry/polygon.py @@ -144,28 +144,11 @@ def polygonToSpeckleMesh( def getZaxisTranslation(layer, boundaryPts, dataStorage): #### check if elevation is applied and layer exists: elevationLayer = getElevationLayer(dataStorage) - polygonWkt = dataStorage.project.crs().toWkt() - polygonProj = ( - QgsCoordinateReferenceSystem.fromWkt(polygonWkt) - .toProj() - .replace(" +type=crs", "") - ) - translationValue = None + if elevationLayer is not None: all_arrays, all_mins, all_maxs, all_na = getRasterArrays(elevationLayer) settings_elevation_layer = get_raster_stats(elevationLayer) - ( - xres, - yres, - originX, - originY, - sizeX, - sizeY, - rasterWkt, - rasterProj, - ) = settings_elevation_layer - allElevations = [] for pt in boundaryPts: # posX, posY = reprojectPt( diff --git a/speckle/converter/geometry/utils.py b/speckle/converter/geometry/utils.py index c43debec..a693256f 100644 --- a/speckle/converter/geometry/utils.py +++ b/speckle/converter/geometry/utils.py @@ -127,7 +127,9 @@ def projectToPolygon( return z -def getPolyPtsSegments(geom: Any, dataStorage: "DataStorage", xform): +def getPolyPtsSegments( + geom: Any, dataStorage: "DataStorage", coef: Union[int, None] = None, xform=None +): vertices = [] vertices3d = [] segmList = [] @@ -145,18 +147,27 @@ def getPolyPtsSegments(geom: Any, dataStorage: "DataStorage", xform): pt_iterator = geom.vertices() # get boundary points and segments - pointListLocal = [] + pointListLocalOuter = [] startLen = len(vertices) for i, pt in enumerate(pt_iterator): if ( - len(pointListLocal) > 0 - and pt.x() == pointListLocal[0].x() - and pt.y() == pointListLocal[0].y() - ): # don't repeat 1st point + len(pointListLocalOuter) > 0 + and pt.x() == pointListLocalOuter[0].x() + and pt.y() == pointListLocalOuter[0].y() + ): + # don't repeat 1st point pass + elif coef is None: + pointListLocalOuter.append(pt) else: - pointListLocal.append(pt) - for i, pt in enumerate(pointListLocal): + if i % coef == 0: + pointListLocalOuter.append(pt) + else: + # don't add points, which are in-between specified step (coeff) + # e.g. if coeff=5, we skip ponts 1,2,3,4, but add points 0 and 5 + pass + + for i, pt in enumerate(pointListLocalOuter): x, y = apply_pt_offsets_rotation_on_send(pt.x(), pt.y(), dataStorage) vertices.append([x, y]) try: @@ -182,18 +193,29 @@ def getPolyPtsSegments(geom: Any, dataStorage: "DataStorage", xform): pt_list = list(pt_iterator) pointListLocal = [] startLen = len(vertices) + for i, pt in enumerate(pt_list): if ( len(pointListLocal) > 0 and pt.x() == pointListLocal[0].x() and pt.y() == pointListLocal[0].y() - ): # don't repeat 1st point + ): + # don't repeat 1st point continue - elif [ - pt.x(), - pt.y(), - ] not in vertices: # in case it's not the inner part of geometry - pointListLocal.append(pt) + elif pt not in pointListLocalOuter: + # make sure it's not already included in the outer part of geometry + + if coef is None or len(pt_list) / coef < 5: + # coef was calculated by the outer ring. + # We need to make sure inner ring will have at least 4 points, otherwise ignore coeff. + pointListLocal.append(pt) + else: + if i % coef == 0: + pointListLocal.append(pt) + else: + # don't add points, which are in-between specified step (coeff) + # e.g. if coeff=5, we skip ponts 1,2,3,4, but add points 0 and 5 + pass if len(pointListLocal) > 2: holes.append( @@ -322,14 +344,14 @@ def to_triangles(data: dict, attempt: int = 0) -> Tuple[Union[dict, None], int]: def triangulatePolygon( - geom: Any, dataStorage: "DataStorage", xform=None + geom: Any, dataStorage: "DataStorage", coef: Union[int, None] = None, xform=None ) -> Tuple[dict, Union[List[List[float]], None], int]: try: # import triangle as tr vertices = [] # only outer segments = [] # including holes holes = [] - pack = getPolyPtsSegments(geom, dataStorage, xform) + pack = getPolyPtsSegments(geom, dataStorage, coef, xform) vertices, vertices3d, segments, holes = pack @@ -432,13 +454,20 @@ def fix_orientation( index = k + 1 if k == len(polyBorder) - 1: index = 0 - pt = polyBorder[k * coef] - pt2 = polyBorder[index * coef] - if isinstance(pt, Point) and isinstance(pt2, Point): - sum_orientation += (pt2.x - pt.x) * (pt2.y + pt.y) # if Speckle Points - else: - sum_orientation += (pt2.x() - pt.x()) * (pt2.y() + pt.y()) # if QGIS Points + try: + pt = polyBorder[k * coef] + pt2 = polyBorder[index * coef] + + if isinstance(pt, Point) and isinstance(pt2, Point): + sum_orientation += (pt2.x - pt.x) * (pt2.y + pt.y) # if Speckle Points + else: + sum_orientation += (pt2.x() - pt.x()) * ( + pt2.y() + pt.y() + ) # if QGIS Points + except IndexError: + break + if positive is True: if sum_orientation < 0: polyBorder.reverse() @@ -735,11 +764,16 @@ def getArcNormal(poly: Arc, midPt: Point, dataStorage) -> Union[Vector, None]: def apply_pt_offsets_rotation_on_send( x: float, y: float, dataStorage -) -> Tuple[Union[float, None]]: # on Send +) -> Tuple[float, float]: # on Send try: offset_x = dataStorage.crs_offset_x offset_y = dataStorage.crs_offset_y rotation = dataStorage.crs_rotation + if offset_x == offset_y == rotation == 0: + return x, y + if offset_x is None and offset_y is None and rotation is None: + return x, y + if offset_x is not None and isinstance(offset_x, float): x -= offset_x if offset_y is not None and isinstance(offset_y, float): @@ -757,7 +791,7 @@ def apply_pt_offsets_rotation_on_send( return x, y except Exception as e: logToUser(e, level=2, func=inspect.stack()[0][3]) - return None, None + raise e def transform_speckle_pt_on_receive(pt_original: Point, dataStorage) -> Point: diff --git a/speckle/converter/layers/symbology.py b/speckle/converter/layers/symbology.py index 58cec088..32c26538 100644 --- a/speckle/converter/layers/symbology.py +++ b/speckle/converter/layers/symbology.py @@ -606,20 +606,16 @@ def rendererToSpeckle( redBand = renderer.redBand() greenBand = renderer.greenBand() blueBand = renderer.blueBand() - redContrast = ( - redMin - ) = ( - redMax - ) = ( - greenContrast - ) = greenMin = greenMax = blueContrast = blueMin = blueMax = None + redContrast = redMin = redMax = greenContrast = greenMin = greenMax = ( + blueContrast + ) = blueMin = blueMax = None try: redContrast = ( renderer.redContrastEnhancement().contrastEnhancementAlgorithm() ) redMin = renderer.redContrastEnhancement().minimumValue() redMax = renderer.redContrastEnhancement().maximumValue() - except: + except: # AttributeError: 'NoneType' object has no attribute 'contrastEnhancementAlgorithm' pass try: greenContrast = ( @@ -627,7 +623,7 @@ def rendererToSpeckle( ) greenMin = renderer.greenContrastEnhancement().minimumValue() greenMax = renderer.greenContrastEnhancement().maximumValue() - except: + except: # AttributeError pass try: blueContrast = ( @@ -635,7 +631,7 @@ def rendererToSpeckle( ) blueMin = renderer.blueContrastEnhancement().minimumValue() blueMax = renderer.blueContrastEnhancement().maximumValue() - except: + except: # AttributeError pass layerRenderer.update( { diff --git a/speckle/converter/layers/utils.py b/speckle/converter/layers/utils.py index e0cba5b9..c8ba49f9 100644 --- a/speckle/converter/layers/utils.py +++ b/speckle/converter/layers/utils.py @@ -517,6 +517,7 @@ def getClosestIndex(x): def getArrayIndicesFromXY(settings, x, y): + """Get cell x&y incices (and remainders) on a given layer, from absolute XY coordinates.""" resX, resY, minX, minY, sizeX, sizeY, wkt, proj = settings index2 = (x - minX) / resX index1 = (y - minY) / resY @@ -545,47 +546,10 @@ def getArrayIndicesFromXY(settings, x, y): return ind1, ind2, remainder1, remainder2 -def getHeightWithRemainderFromArray(height_array, texture_transform, ind1, ind2): - if texture_transform is True: - z = height_array[ind1][ind2] - r""" - # TODO - indexFloor1 = index1_0 if ind1>0 else ind1 - indexFloor2 = index2_0 if ind2>0 else ind2 - - rem1 = remainder1 if ind1>0 else remainder1_0 - rem2 = remainder2 if ind2>0 else remainder2_0 - - rem1_share = rem1/ (rem1 + rem2) - rem2_share = rem2/ (rem1 + rem2) - z = height_array[indexFloor1][indexFloor2] + (height_array[ind1][ind2] - height_array[indexFloor1][ind2])*rem1*rem1_share + (height_array[ind1][ind2] - height_array[ind1][indexFloor2])*rem2*rem2_share - """ - else: - z = height_array[ind1][ind2] - try: - val = float(z) - except: - val = None - return val - - -def getXYofArrayPoint( - settings, indexX, indexY, selectedLayer, elevationLayer, dataStorage -): - resX, resY, minX, minY = settings - x = minX + resX * indexX - y = minY + resY * indexY - # newX, newY = reprojectPt(x, y, wkt, proj, targetWKT, targetPROJ) - reprojected_pt = transform.transform( - dataStorage.project, - QgsPointXY(x, y), - selectedLayer.crs(), - elevationLayer.crs(), - ) - newX = reprojected_pt.x() - newY = reprojected_pt.y() - - return newX, newY +def getXYofArrayPoint(rasterResXY, minX, minY, indexX, indexY): + x = minX + rasterResXY[0] * indexX + y = minY + rasterResXY[1] * indexY + return x, y def isAppliedLayerTransformByKeywords( @@ -631,6 +595,8 @@ def isAppliedLayerTransformByKeywords( def getElevationLayer(dataStorage): elevationLayer = dataStorage.elevationLayer + if elevationLayer is None: + return None try: # check if layer was not deleted name = elevationLayer.name() @@ -657,7 +623,7 @@ def get_raster_stats(rasterLayer): sizeX, sizeY = (band.ReadAsArray().shape[1], band.ReadAsArray().shape[0]) return xres, yres, originX, originY, sizeX, sizeY, rasterWkt, rasterProj - except: + except Exception as e: return None, None, None, None, None, None, None, None @@ -666,17 +632,6 @@ def getRasterArrays(elevationLayer): try: elevationSource = gdal.Open(elevationLayer.source(), gdal.GA_ReadOnly) - settings_elevation_layer = get_raster_stats(elevationLayer) - ( - xres, - yres, - originX, - originY, - sizeX, - sizeY, - rasterWkt, - rasterProj, - ) = settings_elevation_layer all_arrays = [] all_mins = [] diff --git a/speckle/utils/project_vars.py b/speckle/utils/project_vars.py index 6da45b6c..d3fabbe3 100644 --- a/speckle/utils/project_vars.py +++ b/speckle/utils/project_vars.py @@ -283,10 +283,11 @@ def setProjectReferenceSystem(dataStorage: DataStorage, dockwidget=None): try: from qgis.core import QgsCoordinateReferenceSystem + r""" newCrsString = ( "+proj=tmerc +ellps=WGS84 +datum=WGS84 +units=m +no_defs +lon_0=" + str(dataStorage.custom_lon) - + " lat_0=" + + " +lat_0=" + str(dataStorage.custom_lat) + " +x_0=0 +y_0=0 +k_0=1" ) @@ -294,8 +295,10 @@ def setProjectReferenceSystem(dataStorage: DataStorage, dockwidget=None): newCrsString ) # fromWkt(newProjWkt) validate = QgsCoordinateReferenceSystem().createFromProj(newCrsString) - wkt = newCrs.toWkt() + """ + wkt = f'PROJCS["SpeckleCRS_latlon_{dataStorage.custom_lat}_{dataStorage.custom_lon}", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["Degree", 0.0174532925199433]], PROJECTION["Transverse_Mercator"], PARAMETER["False_Easting", 0.0], PARAMETER["False_Northing", 0.0], PARAMETER["Central_Meridian", {dataStorage.custom_lon}], PARAMETER["Scale_Factor", 1.0], PARAMETER["Latitude_Of_Origin", {dataStorage.custom_lat}], UNIT["Meter", 1.0]]' + validate = QgsCoordinateReferenceSystem().createFromWkt(wkt) newCRSfromWkt = QgsCoordinateReferenceSystem.fromWkt(wkt) if validate: diff --git a/speckle_qgis.py b/speckle_qgis.py index 7531351a..9693fdf4 100644 --- a/speckle_qgis.py +++ b/speckle_qgis.py @@ -496,7 +496,8 @@ def onSend(self, message: str): units = str(QgsUnitTypes.encodeUnit(projectCRS.mapUnits())) self.dataStorage.latestActionUnits = units try: - units = get_units_from_string(units) + units_class = get_units_from_string(units) + units = units_class.value except SpeckleInvalidUnitException: units = "none" self.dataStorage.currentUnits = units diff --git a/tests/unit/converter/layer_tests/test_layer_utils.py b/tests/unit/converter/layer_tests/test_layer_utils.py index 0395c0ce..2d060a8b 100644 --- a/tests/unit/converter/layer_tests/test_layer_utils.py +++ b/tests/unit/converter/layer_tests/test_layer_utils.py @@ -11,7 +11,6 @@ reprojectPt, getClosestIndex, getArrayIndicesFromXY, - getHeightWithRemainderFromArray, getXYofArrayPoint, isAppliedLayerTransformByKeywords, getElevationLayer, diff --git a/tests/unit/utils/test_validation.py b/tests/unit/utils/test_validation.py index fc9d5bbc..015c7b80 100644 --- a/tests/unit/utils/test_validation.py +++ b/tests/unit/utils/test_validation.py @@ -116,8 +116,11 @@ def test_validateBranch(stream): assert isinstance(result, Branch) -def test_validateBranch_no_commits(stream_fe1): - branch_name = "main" +def test_validateBranch_no_commits(): + sample_wrapper = StreamWrapper("https://latest.speckle.dev/streams/7117052f4e") + sample_wrapper.get_client() + stream_fe1 = sample_wrapper._client.stream.get(id=sample_wrapper.stream_id) + branch_name = "empty_branch" result = validateBranch(stream_fe1, branch_name, checkCommits=True, dockwidget=None) assert result is None