Add Geoid to DEM-Export

[insanejules]

I am trying to give my DEM a projection and a geoid. When I use the same geoid file in Metashape Pro it works. However, when I try to use it in the Python API I get the following error:

Code: [Select]

proj = Metashape.OrthoProjection()
proj.type = Metashape.OrthoProjection.Type.Planar
#Metashape.CoordinateSystem.addGeoid("de_bkg_GCG2016v2023.tif")
proj.crs = Metashape.CoordinateSystem("EPSG:25832")
chunk.buildDem(source_data=Metashape.PointCloudData, interpolation = Metashape.Interpolation.DisabledInterpolation)
doc.save()

if chunk.elevation:
    chunk.exportRaster('output/dem.tif', source_data= Metashape.ElevationData, projection = proj)
ERROR:

Code: [Select]

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[31], line 4
      2 proj.type = Metashape.OrthoProjection.Type.Planar
      3 #Metashape.CoordinateSystem.addGeoid("de_bkg_GCG2016v2023.tif")
----> 4 proj.crs = Metashape.CoordinateSystem("EPSG:25832")
      5 chunk.buildDem(source_data=Metashape.PointCloudData, interpolation = Metashape.Interpolation.DisabledInterpolation)
      6 doc.save()

AttributeError: 'Metashape.Metashape.OrthoProjection' object attribute 'crs' is read-only

[Alexey Pasumansky]

Hello insanejules,

Do you observe that AttributeError, if you re-start Metashape instance and try to use same code?
I assume that there could be some incorrect assignment in the same instance of Metashape that broke the interpreter. So at first you need to solve the problem, with this code.

Then, if you need to use compound coordinate system for DEM generation and/or export, you should assign it to proj.crs - either use available systems from EPSG registry (like ETRS89 / UTM zone 32N + DHHN92 height - EPSG::5555) or to use wkt line for assignment, if there the desired coordinate system is not available by EPSG code. For example:

Code: [Select]

wkt = 'COMPD_CS["ETRS89 / UTM zone 32N + DHHN2016 height",
  PROJCS["ETRS89 / UTM zone 32N",
  GEOGCS["ETRS89",
    DATUM["European Terrestrial Reference System 1989 ensemble",
      SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
      TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]], AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator",AUTHORITY["EPSG","9807"]],
    PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","25832"]],
  VERT_CS["DHHN2016 height",
    VERT_DATUM["Deutsches Haupthoehennetz 2016",2005,AUTHORITY["EPSG","1170"]],
    UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","7837"]]]'
proj.crs = Metashape.CoordinateSystem(wkt)

Since you say that you have the procedure working via GUI, then you can save the compound coordinate system used to PRJ file, open it as text file and copy the definition to your script as a string.

Also you would need to uncomment Metashape.CoordinateSystem.addGeoid("de_bkg_GCG2016v2023.tif") line in the beginning of the script (prior any coordinate system assignments are used) and pass the path to the corresponding geoid file to the function.

[insanejules]

Thank you very much Alexey, that wkt-definition works great !

Here is the code for those who have the same problem:

Code: [Select]

proj = Metashape.OrthoProjection()
Metashape.CoordinateSystem.addGeoid("de_bkg_GCG2016v2023.tif")

wkt = 'COMPD_CS["ETRS89 / UTM zone 32N + DHHN2016 height", PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European Terrestrial Reference System 1989 ensemble",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9102"]], AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","25832"]],VERT_CS["DHHN2016 height",VERT_DATUM["Deutsches Haupthoehennetz 2016",2005,AUTHORITY["EPSG","1170"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","7837"]]]'


proj.crs = Metashape.CoordinateSystem(wkt)
chunk.buildDem(source_data=Metashape.PointCloudData, interpolation = Metashape.Interpolation.DisabledInterpolation, projection = proj)
doc.save()
chunk.crs = Metashape.CoordinateSystem(wkt)
print(chunk.crs)
if chunk.elevation:
    chunk.exportRaster('output/dem.tif', source_data= Metashape.ElevationData, projection = proj, nodata_value=999)