Skip to content

geoid issue-19 #20

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

Jump to bottom
Open
wants to merge 27 commits into
base: main
Choose a base branch
from
Open

geoid issue-19 #20

wants to merge 27 commits into from

Conversation

brentwilder
Copy link
Member

@brentwilder brentwilder commented Jan 21, 2025 •
edited
Loading

Fixes hard coding of CRS for geoid conversion. (Closes #19 )

  • CRS of correct projection is pulled directly from the input laz

@brentwilder brentwilder self-assigned this Jan 21, 2025
@jomey
Copy link
Member

jomey commented Jan 21, 2025

Looking at the initial issue and the implementation of this PR, I am going back and forth on the best approach. On the one side it is nice to have some sort of automation for converting between CRS and the other is that we are assuming that the point cloud is the preferred/correct one.

You mentioned pulling this logic out of the pipeline and have it as a separate tool/script, which I am currently leaning towards.

For this PR, I would suggest that we compare the reference DEM and the point cloud and raise an error in the workflow when there is a difference. Then we could give the user the exact conversion call for GDAL if they opt to convert the reference DEM. Basically print the call on line 74

Making the user aware of this detail is important and then have ice-roads as the tool to assist the user in getting everything done.

Thoughts?

@jomey jomey requested a review from naniciafone January 21, 2025 23:38
@shadoneel
Copy link
Collaborator

shadoneel commented Jan 22, 2025 via email

@brentwilder
Copy link
Member Author

@shadoneel But we do ask for this in the README (see below),

-g geoid         Is the reference DEM CRS orthometric (geoid height)? Will be auto set to True if you don't supply a DEM [Default: False]

Users input either a True or False at the start.

imo a GUI may be a little over-engineered, and would not be ideal for running on Borah or other linux cluster.

@jomey My main issue with pulling it out, is that it's not just GDAL , since GDAL does not have a geoid conversion. We would be asking users to compute it using Ames Stereo. Which could be less accessible. It's possible that keeping this in as is, would aid non-technical users in easily performing the geoid conversion?

Alternatively, we can do what I was suggesting and just remove the -g flag. We can have an exception that looks for if reference-DEM is in geoid. And we can refer them to run another script that we make. But, this whole process could be longer, and potentially not worth the headache if we can just ask the user to be aware of what vertical datum their reference-DEM is in?

@jomey
Copy link
Member

jomey commented Jan 22, 2025

I am still voting for referring the user to a separate script to properly prepare the reference DEM. Also thinking from a repeatability perspective, where for instance in our case we would convert the reference DEM every time we process a point cloud. This adds unnecessary time each time a workflow is executed. It already takes quite a while for one run and saving time should be desirable.

@shadoneel
Copy link
Collaborator

shadoneel commented Jan 22, 2025 via email

@brentwilder
Copy link
Member Author

@shadoneel @jomey Putting these two ideas together, I would argue we need to have a separate utility script even more so. This would allow for more flexibility (perhaps allowing choice of a GDAL or Ames method), adding in if its in feet vs meter, allowing for users to assess the output stats to verify, repeatability so we are not computing the same thing over again?

@brentwilder
Copy link
Member Author

And I can make it clear in the documentation/README, that this geoid-conversion-tool should be used first.

We can also add another Exception/Error, to catch this early on in the main pipeline.

@jomey
Copy link
Member

jomey commented Jan 22, 2025

Tagging @naniciafone to chime in how she processed the 3DEP data. Think there is also a PDAL+GDAL way that uses the 3DEP point clouds.

@naniciafone
Copy link

PDAL was able to perform the transformation from NAVD 1988 GEOID 12B to WGS ellipsoid with the reprojection filter... I thought I was going to need to download a gtx file but I didn't need to.

@brentwilder
Copy link
Member Author

brentwilder commented Jan 22, 2025 •
edited
Loading

Good to know about the ease of the PDAL method for point clouds @naniciafone and @jomey ... Which could also be added to this geoid-helper script?

Are we in consensus then to take out the -g flag and create a sort of standalone tool that can help guide some of these approaches?

EDIT: with the idea being similar to Shad's comments that we are focusing on User input, and trying to simply bring down the technical barrier on this part?

@jomey
Copy link
Member

jomey commented Jan 22, 2025 •
edited
Loading

Good to know about the ease of the PDAL method for point clouds @naniciafone and @jomey ... Which could also be added to this geoid-helper script?

I would like to see Nani's workflow at least in the example. Think this is very useful to have as a reference if a user decides to go the point cloud route.

Are we in consensus then to take out the -g flag and create a sort of standalone tool that can help guide some of these approaches?

+1 from my side

EDIT: with the idea being similar to Shad's comments that we are focusing on User input, and trying to simply bring down the technical barrier on this part?

jomey
jomey reviewed Jan 24, 2025
View reviewed changes
jomey
jomey reviewed Jan 24, 2025
View reviewed changes
jomey
jomey reviewed Jan 24, 2025
View reviewed changes
jomey
jomey reviewed Jan 24, 2025
View reviewed changes
@brentwilder
Copy link
Member Author

Should need a little more testing perhaps? but i have exhausted all of the data I have on my computer (flipping back and forth between WGS and NAD). But I believe I have done it? It seems to be fully dummy proof and does not allow you to do any sort of geoid transformations that are impossible.

It also takes in either a raster (using ASP), or

it can take in a point cloud (using PDAL). I leverage the fact we have the navd88.tif already local due to the ASP install. And so this is called in the PDAL pipeline to ensure consistency when performing the transformation.

@jomey
Copy link
Member

jomey commented Jan 24, 2025

Should need a little more testing perhaps?

@naniciafone - Want to give this a spin with your current use case in Kamas?

@brentwilder
Copy link
Member Author

On the one side it is nice to have some sort of automation for converting between CRS and the other is that we are assuming that the point cloud is the preferred/correct one.

Also, just to follow up on your earlier point @jomey , we have taken out this assumption too. And so now, we just ask the user to provide the EPSG code (e.g 32611) that they wish their data to be in. This should hopefully remove the ambiguity in which is the preferred/correct one.

brentwilder
brentwilder commented Jan 25, 2025
View reviewed changes
brentwilder
brentwilder commented Jan 25, 2025
View reviewed changes
@brentwilder
Copy link
Member Author

@naniciafone perhaps the easiest and most flexible solution for us would be to have another argument . like an "override" argument that mostly would just be used for USGS 3DEP lidar. that could perform geoid->ellipsoid by ignoring the datum? ... Seems a little dangerous just because it expects the user to know exactly. But i think this is required due to the ambiguity?

@brentwilder
Copy link
Member Author

Okay, I have pushed a new version if you wanted to test! Updates:

  • You may put a -u flag to override the vertical datum info in the laz file. No argument required.
  • download_dem() has been removed. User DEM is now required.
  • per @Surfix suggestion the name is now changed too to be , transform_data.py.

Example usage would be:

python ./scripts/transform_data.py test.laz -e 32611 -t to_ellipsoid -a /asp/example -u

@naniciafone
Copy link

yeah I can try it again! Question for my own understanding- does the tool automatically convert to WGS84 (versus NAD83/ another ellipsoid)? And does it assume a geoid model? If the geoid model is not specified in the metadata, I assume I'd need to specify. I've only done transformations with PDAL/PROJ so not sure how ASP handles it.

@brentwilder
Copy link
Member Author

brentwilder commented Mar 21, 2025 •
edited
Loading

Yeah that’s a great question. And right now it’s super simple, and just assumes ‘navd88.tif’. Which is supplied in ASP. We use this same model for both point cloud and raster data.

I’m realizing as I am writing this , that I need to check if we need to use the WGS model in here too (EGM2008).

Update: this is a great point Nani I need to think more about this.

@brentwilder
Copy link
Member Author

brentwilder commented Mar 21, 2025 •
edited
Loading

Yeah so I am missing cases here and the code is currently flawed.. It assumes user wants to switch from either NAD83 to WGS or visa versa.

But you can stay within one of those datums and still change geoid to ellipsoid.

I will continue to iterate and figure out the best way to do this but it will take some time.

@naniciafone
Copy link

naniciafone commented Mar 21, 2025 •
edited
Loading

If it's helpful, this is what I ended up doing... with much hand holding from Abhi @ NCALM, who visited for the SNOWWI flights :) :

  1. I ended up downloading the AWS data in EPSG 6340 (NAD83/ UTM 11N). So 6340 height, NAVD88 (EPSG 5703)heights. I was told that you have to convert NAVD88 to NAD83 BEFORE converting to WGS84. NAVD88 + Geoid height = NAD83.
  2. Used a JSON to convert NAD83/UTM to geographic coordinates AND perform that calculation to get from NAVD88 to NAD83. I initially used a .gtx grid file but PROJ wasn't happy with that, so I downloaded a .tif from here instead- https://cdn.proj.org/ :
[
	{
	"type" : "readers.las",
	"filename" : "pc_sample.laz"
	},
	{
	"type":"filters.projpipeline",
	"coord_op":"+proj=pipeline +step +inv +proj=utm +zone=11 +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"
	},
	{
	"type":"filters.projpipeline",	
	"coord_op":"+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_noaa_geoid09_conus.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"
	},
	{
	"type":"filters.projpipeline",	
	"coord_op":"+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=11 +ellps=GRS80"
	},
	{
	"type" : "writers.las",
	"filename" : "pc_ellipsoid.laz",
	"forward": "all",
	"minor_version":"2",
	"a_srs": "EPSG:6340",
	"system_id":"RIEGL-580-II",
	"software_id":"TERRASCAN_PDAL"
	}
]
  1. Used a second json to convert from NAD83 to WGS84 AND WGS84 geographic coordinated --> WGS84/UTM11N:
[
	{
	"type" : "readers.las",
	"filename" : "pc_ellipsoid.laz"
	},
	{
	"type": "filters.projpipeline",
	"coord_op": "  +proj=pipeline +step +inv +proj=utm +zone=11 +ellps=GRS80 +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138 +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006 +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05 +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1"
	},
	{
	"type": "filters.projpipeline",
	"coord_op":"+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=11 +ellps=WGS84"
	},
	{
	"type" : "writers.las",
	"filename" : "pc_WGS84.laz",
	"forward": "all",
	"minor_version":"2",
	"a_srs": "EPSG:32611",
	"system_id":"RIEGL-580-II"
	}
]

Probably could've done one json but it was helpful to look at my point samples along the way. The 3DEP and RaLiPod were 15m off to start. After converting to NAD83 with the first json, the points were about 30cm offset. The second json removed that last 30cm offset.

@naniciafone
Copy link

naniciafone commented Mar 21, 2025 •
edited
Loading

PDAL was able to perform the transformation from NAVD 1988 GEOID 12B to WGS ellipsoid with the reprojection filter... I thought I was going to need to download a gtx file but I didn't need to.

I haven't been able to replicate this since I moved to the AWS downloaded point cloud, which as discussed does not have encoded vertical CRS data... when I posted this I was using tiled data downloaded from the National Map, which was a headache for a lot of reasons. That data did have encoded vertical CRS info, however.... I've also learned a ton about transformations since I posted this and want to verify that I actually got the transformation I thought I did.

@brentwilder
Copy link
Member Author

brentwilder commented Mar 23, 2025 •
edited
Loading

So a few things as I am looking at these examples (thank you for providing).

  • To convert between ellipsoids (NAD83 <-> WGS), you must include 7-parameter Helmert transformation parameters which are space and time dependent, correct? Too many possibilities to make like a standard tool.
  • Further, I'm now having a hard time thinking how we wrap this into an actual tool>? There are just so many variations that someone may need for potential transformation....
    • Is it possible that we may just point to resources ?? Like a Folder in our repo with this json example for point clouds using PDAL? And then I could craft a short example using ASP's dem_geoid (https://stereopipeline.readthedocs.io/en/latest/tools/dem_geoid.html) for rasters? And then in the readme we can just briefly state that Reference DEM and Point Clouds must be in same horizontal/vertical reference system as we already have?

EDIT: I have gone ahead and updated it to be more of a "library" of sorts that users may start from, versus a formal tool. I honestly think this will be more helpful in the long-term and perhaps @naniciafone like as you use them you could put them here for others (see contrib/datum_transformation_examples/)? Let me know what you think.

@jomey
Copy link
Member

jomey commented Mar 24, 2025

I like the idea of an example folder instead of the tool (for the reasons Brent mentioned). We could publish the PDAL JSON pipelines that handle 3DEP and any other ones that you crated, Brent.
Then the only logic remaining for Ice-Roads would be to check the information before starting and if there is not a match, error out and link to the example folder as a first help to users?

@jomey
Copy link
Member

jomey commented Mar 24, 2025

So a few things as I am looking at these examples (thank you for providing).

* To convert between ellipsoids (NAD83 <-> WGS), you must include 7-parameter Helmert transformation parameters which are space and time dependent, correct? Too many possibilities to make like a standard tool.

Going back to this question:
https://www.ngs.noaa.gov/CORS/news/historical_helmert.shtml

@naniciafone
Copy link

Happy to add my own examples to that folder/ contribute to the "resources" folder. I think that's a good way to go given the nuances of geoid-ellipsoid transformations.

@shadoneel
Copy link
Collaborator

shadoneel commented Mar 24, 2025 via email

@naniciafone
Copy link

naniciafone commented Mar 24, 2025 •
edited
Loading

I think the way we could make it a tool is by helping the user build the json with all PROJ strings needed for transformations. I got the proj strings through, for example, projinfo -s "EPSG:5703" -t "EPSG:6340" -o PROJ -q. The thing is, the user still needs to know all steps for the transformation, because I think PROJ will only output a single string... I'm not totally sure why that is.

@brentwilder
Copy link
Member Author

@naniciafone I agree and I came to a similar conclusion when I was looking at that too.

That's why I think building out an example library would be most helpful. And especially for Mores crew, and any of your more recent projects @naniciafone and @shadoneel to account for the largest use cases!

If you want to email the jsons to mee I can just go ahead and add them to the https://github.com/cryogars/ice-road-copters/tree/issue-19-geoid/contrib/datum_transformation_examples , and do another push to this PR.

Let me know!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@jomey jomey jomey left review comments

@shadoneel shadoneel Awaiting requested review from shadoneel

@ZachHoppinen ZachHoppinen Awaiting requested review from ZachHoppinen

@naniciafone naniciafone Awaiting requested review from naniciafone

At least 1 approving review is required to merge this pull request.

Assignees

@brentwilder brentwilder

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

download_dem() via py3dep (USGS 3DEP) incorrectly assumes ellipsoid CRS hard coded in Geoid conversion
5 participants
@brentwilder @jomey @shadoneel @naniciafone @Surfix