Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improved Geo VLR Metadata #11

Merged
merged 2 commits into from
Sep 24, 2024
Merged

Improved Geo VLR Metadata #11

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
575b946
Add WKT vlr
tayden Sep 24, 2024
e0ee636
Improve Geo VLR attributes
tayden Sep 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 34 additions & 5 deletions las_trx/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
from datetime import date
from typing import Optional

from csrspy.enums import CoordType, Reference, VerticalDatum
from csrspy.utils import date_to_decimal_year
from pydantic import BaseModel
from pyproj.crs import (
CRS,
Expand All @@ -14,9 +16,6 @@
from pyproj.crs.coordinate_operation import UTMConversion
from pyproj.crs.coordinate_system import Cartesian2DCS

from csrspy.enums import CoordType, Reference, VerticalDatum
from csrspy.utils import date_to_decimal_year


class CSRSPYConfig(BaseModel):
s_ref_frame: Reference
Expand All @@ -38,9 +37,39 @@ class TrxVd(str, enum.Enum):

@property
def vertical_crs(self) -> Optional[VerticalCRS]:
if self in [TrxVd.CGG2013A, TrxVd.CGG2013]:
# Seems to be no distinction between 2013 and 2013a
if self == TrxVd.CGG2013:
return VerticalCRS.from_epsg(6647)
elif self == TrxVd.CGG2013A:
return VerticalCRS.from_dict(
{
"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json",
"type": "VerticalCRS",
"name": "CGVD2013(CGG2013a) height",
"datum": {
"type": "VerticalReferenceFrame",
"name": "Canadian Geodetic Vertical Datum of 2013 (CGG2013a)",
},
"coordinate_system": {
"subtype": "vertical",
"axis": [
{
"name": "Gravity-related height",
"abbreviation": "H",
"direction": "up",
"unit": "metre",
}
],
},
"scope": "Geodesy, engineering survey, topographic mapping.",
"area": "Canada - onshore and offshore - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon.",
"bbox": {
"south_latitude": 38.21,
"west_longitude": -141.01,
"north_latitude": 86.46,
"east_longitude": -40.73,
},
}
)
elif self == TrxVd.HT2_2010v70:
return VerticalCRS.from_epsg(5713)
return None
Expand Down
153 changes: 11 additions & 142 deletions las_trx/vlr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,65 +1,24 @@
"""
Copied and modified from laspy known vlr module.

Provides easier and corrected writing of geo tags
Provides easier and corrected writing of geotags

Should switch to the official laspy version if they get this module cleaned up a bit.
"""

import ctypes
from typing import Type, TypeVar
from typing import Type

import pyproj
from laspy.vlrs import BaseKnownVLR
from laspy.vlrs.known import (
GeoAsciiParamsVlr,
GeoKeyDirectoryVlr,
GeoKeyEntryStruct,
GeoAsciiParamsType,
GeoKeyDirectoryType,
)

NULL_BYTE = b"\0"

GeoKeyDirectoryType = TypeVar("GeoKeyDirectoryType", bound="GeoKeyDirectoryVlr")
GeoAsciiParamsType = TypeVar("GeoAsciiParamsType", bound="GeoAsciiParamsVlr")


class GeoKeyEntryStruct(ctypes.LittleEndianStructure):
_pack_ = 1
_fields_ = [
("id", ctypes.c_uint16),
("tiff_tag_location", ctypes.c_uint16),
("count", ctypes.c_uint16),
("value_offset", ctypes.c_uint16),
]

def __init__(self, id=0, tiff_tag_location=0, count=0, value_offset=0):
super().__init__(
id=id,
tiff_tag_location=tiff_tag_location,
count=count,
value_offset=value_offset,
)

@staticmethod
def size():
return ctypes.sizeof(GeoKeysHeaderStructs)

def __repr__(self):
return "<GeoKey(Id: {}, Location: {}, count: {}, offset: {})>".format(
self.id, self.tiff_tag_location, self.count, self.value_offset
)


class GeoAsciiParamsVlr(BaseKnownVLR):
def __init__(self):
super().__init__(description="GeoTIFF GeoAsciiParamsTag")
self.strings = []

def parse_record_data(self, record_data):
st = [s.decode("ascii") for s in record_data.split(NULL_BYTE)]
if len(st) > 1:
self.strings = st[0][:-1].split("|")
else:
self.strings = []

def record_data_bytes(self):
return NULL_BYTE.join([self.ascii_params, b""])

class TrxGeoAsciiParamsVlr(GeoAsciiParamsVlr):
def record_from_crs(self, crs):
if crs.geodetic_crs is not None:
self.strings.append(crs.geodetic_crs.name)
Expand All @@ -74,93 +33,14 @@ def record_from_crs(self, crs):
_, vd_name = crs.name.split(" + ")
self.strings.append(vd_name)

@property
def ascii_params(self):
return "".join([s + "|" for s in self.strings]).encode("ascii")

@classmethod
def from_crs(cls: Type[GeoAsciiParamsType], crs: pyproj.CRS) -> GeoAsciiParamsType:
self = cls()
self.record_from_crs(crs)
return self

def __repr__(self):
return "<GeoAsciiParamsVlr({})>".format(self.strings)

@staticmethod
def official_user_id():
return "LASF_Projection"

@staticmethod
def official_record_ids():
return (34737,)


class GeoKeysHeaderStructs(ctypes.LittleEndianStructure):
_pack_ = 1
_fields_ = [
("key_directory_version", ctypes.c_uint16),
("key_revision", ctypes.c_uint16),
("minor_revision", ctypes.c_uint16),
("number_of_keys", ctypes.c_uint16),
]

def __init__(
self,
key_directory_version=1,
key_revision=1,
minor_revision=0,
number_of_keys=0,
):
super().__init__(
key_directory_version=key_directory_version,
key_revision=key_revision,
minor_revision=minor_revision,
number_of_keys=number_of_keys,
)

@staticmethod
def size():
return ctypes.sizeof(GeoKeysHeaderStructs)

def __repr__(self):
return "<GeoKeysHeader(vers: {}, rev:{}, minor: {}, num_keys: {})>".format(
self.key_directory_version,
self.key_revision,
self.minor_revision,
self.number_of_keys,
)


class GeoKeyDirectoryVlr(BaseKnownVLR):
def __init__(self):
super().__init__(description="GeoTIFF GeoKeyDirectoryTag")
self.geo_keys_header = GeoKeysHeaderStructs()
self.geo_keys = [GeoKeyEntryStruct()]

def parse_record_data(self, record_data):
record_data = bytearray(record_data)
header_data = record_data[: ctypes.sizeof(GeoKeysHeaderStructs)]
self.geo_keys_header = GeoKeysHeaderStructs.from_buffer(header_data)
self.geo_keys = []
keys_data = record_data[GeoKeysHeaderStructs.size() :]
num_keys = (
len(record_data[GeoKeysHeaderStructs.size() :]) // GeoKeyEntryStruct.size()
)
if num_keys != self.geo_keys_header.number_of_keys:
self.geo_keys_header.number_of_keys = num_keys

for i in range(self.geo_keys_header.number_of_keys):
data = keys_data[
(i * GeoKeyEntryStruct.size()) : (i + 1) * GeoKeyEntryStruct.size()
]
self.geo_keys.append(GeoKeyEntryStruct.from_buffer(data))

def record_data_bytes(self):
b = bytes(self.geo_keys_header)
b += b"".join(map(bytes, self.geo_keys))
return b

class TrxGeoKeyDirectoryVlr(GeoKeyDirectoryVlr):
def record_from_crs(self, crs: pyproj.CRS):
all_keys = {
"GTModelTypeGeoKey": GeoKeyEntryStruct(id=1024, count=1),
Expand Down Expand Up @@ -270,14 +150,3 @@ def from_crs(
self = cls()
self.record_from_crs(crs)
return self

def __repr__(self):
return "<{}({} geo_keys)>".format(self.__class__.__name__, len(self.geo_keys))

@staticmethod
def official_user_id():
return "LASF_Projection"

@staticmethod
def official_record_ids():
return (34735,)
8 changes: 5 additions & 3 deletions las_trx/worker.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,13 @@
import numpy as np
from PyQt6.QtCore import QThread, pyqtSignal as Signal
from laspy import LasHeader
from laspy.vlrs.known import WktCoordinateSystemVlr
from pyproj import CRS

from csrspy import CSRSTransformer
from las_trx.config import TransformConfig
from las_trx.logger import logger
from las_trx.vlr import GeoAsciiParamsVlr, GeoKeyDirectoryVlr
from las_trx.vlr import TrxGeoAsciiParamsVlr, TrxGeoKeyDirectoryVlr

CHUNK_SIZE = 10_000

Expand Down Expand Up @@ -203,8 +204,9 @@ def clear_header_geokeys(header: "LasHeader") -> "LasHeader":


def write_header_geokeys_from_crs(header: "LasHeader", crs: "CRS") -> "LasHeader":
header.vlrs.append(GeoAsciiParamsVlr.from_crs(crs))
header.vlrs.append(GeoKeyDirectoryVlr.from_crs(crs))
header.vlrs.append(TrxGeoAsciiParamsVlr.from_crs(crs))
header.vlrs.append(TrxGeoKeyDirectoryVlr.from_crs(crs))
header.vlrs.append(WktCoordinateSystemVlr(crs.to_wkt()))
logger.debug(f"{header.vlrs=}")
return header

Expand Down