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.

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

First support for ICESat GLAH06 #17

Merged
merged 10 commits into from
Sep 7, 2022
Merged

Conversation

alex-s-gardner
Copy link
Collaborator

@evetion

Here is the first crack at a GLAH06 reader. Two things you might be able to help with:

  1. What is your approach to making nice ascii tables for the function help. I found formatting by hand tricky and didn't manage to get it all

  2. ICESat elevations need to be converted from the TOPEX/POSEIDON to WGS84 the ellipsoid. Do you know of a good way to do this in Julia?

@evetion
Copy link
Owner

evetion commented Aug 30, 2022

That was quick! Nicely done. With regards to your points:

  1. I use https://www.tablesgenerator.com/markdown_tables, you can paste in anything from Excel or similar in there and then export markdown. Any values in the markdown with _ (such as field definitions), need to be escaped with \.
  2. This is a bit of a pain to get right, but I use PROJ to transform from ellipsoid to geoid. There's a function for that (egm2008!). PROJ can also be used to transform between ellipsoids, but this example should be tested and validated, it's very easy to make mistakes here.

Generate a pipeline with PROJ:
❯ projinfo -s "+proj=latlon +a=6378136.299999999813735 +b=6356751.600562936626375 +rf=298.257 +e=0.08181922146 +towgs84=0,0,0 +type=crs" -t EPSG:4979 --3d

Then use it in Julia:

julia> pipe = Proj.proj_create("+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +inv +proj=longlat +a=6378136.3 +rf=298.257 +e=0.08181922146 +step +proj=cart +a=6378136.3 +rf=298.257 +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1")
Ptr{Nothing} @0x00000001254061d0

julia> Proj.proj_trans(pipe, Proj.PJ_FWD, (0,0,0))
4-element Proj.Coord:
  0.0
  0.0
 -0.7000000001862645  # <-- so 70cm lower it seems
 Inf

convert from TOPEX/POSEIDON to WGS84 ellipsoid using Proj.jl
@alex-s-gardner
Copy link
Collaborator Author

@evetion

Once you look this over it should be good to merge

Copy link
Owner

@evetion evetion left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you also write an integration test? See test/runtest.jl for the current suite. Would be good to test that the points function works and returns some actual rows, you can copy the GLAH14 test.

To make the reference file available, add a line like

download_artifact(v"0.1", "GLAH06_634_2131_002_0084_4_01_0001.H5")

src/ICESat/GLAH06.jl Outdated Show resolved Hide resolved
src/ICESat/GLAH06.jl Outdated Show resolved Hide resolved
src/ICESat/GLAH06.jl Outdated Show resolved Hide resolved
pipe = Proj.proj_create("+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +inv +proj=longlat +a=6378136.3 +rf=298.257 +e=0.08181922146 +step +proj=cart +a=6378136.3 +rf=298.257 +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1")

# the values passed to proj_trans() respect the axis order and axis unit of the official definition ( so for example, for EPSG:4326, with latitude first and longitude next, in degrees)
_, _, height_ref = Proj.proj_trans(pipe, Proj.PJ_FWD, (latitude,longitude, height_ref))
Copy link
Owner

@evetion evetion Sep 1, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't work on arrays right? You might need a broadcast and a zip here like: Proj.proj_trans.(pipe, Proj.PJ_FWD, zip(latitude,longitude, height_ref))

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you also write an integration test? See test/runtest.jl for the current suite. Would be good to test that the points function works and returns some actual rows, you can copy the GLAH14 test.

Now implemented

@alex-s-gardner
Copy link
Collaborator Author

all requested fixes and changes implemented

@evetion
Copy link
Owner

evetion commented Sep 5, 2022

I've fixed some small things and removed the following lines

latitude = [x[1] for x in pts]::Vector{Float64}
latitude = [x[1] for x in pts]::Vector{Float64}

as they should be unchanged from their input (only vertical transformation) and allocate extra memory. For the testfile, it seems all shots have quality=false. We might be just unlucky with the granule, but it's something you might check with other granules.

Otherwise I'm happy to merge this.

@alex-s-gardner
Copy link
Collaborator Author

Thanks for the fixes...

I think there are 1.0e-06 differences in latitude because of the slight differences in ellipsoid so it might be best to keep, I'll add back in.

latitude = [x[1] for x in pts]::Vector{Float64}

Let be check on the quality masking

@alex-s-gardner
Copy link
Collaborator Author

Broadcast was missing in quality flag:
(i_numPk .== 1) .&

@evetion
Copy link
Owner

evetion commented Sep 5, 2022

Good catch in the quality broadcast, missed that completely.

Do you expect differences in the latitude for an ellipsoid transformation? What we pass in are float32s, so they're only accurate up to 7 decimals which would be errors of centimeters or less.

@alex-s-gardner
Copy link
Collaborator Author

🤞 we're good to merge

@alex-s-gardner
Copy link
Collaborator Author

Yes, I expect very small differences in the latitude for an ellipsoid transformation (at least I think I have that right).. longitude is unchanged.

But I don't think it matters at all for ICESat data so I'll remove that again

@alex-s-gardner
Copy link
Collaborator Author

latitude = [x[1] for x in pts]::Vector{Float64} is now commented out

@alex-s-gardner
Copy link
Collaborator Author

@evetion you should be good to merge.. then I can start on other items

@evetion evetion merged commit dcf8811 into evetion:master Sep 7, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants