ioos / APIRUS

API for Regular, Unstructured and Staggered model output (or API R US)
Creative Commons Zero v1.0 Universal
2 stars 1 forks source link

Isosurfaces #6

Closed ocefpaf closed 8 years ago

ocefpaf commented 8 years ago

Isosurfaces is one of the methods raised here.

I have a very inefficient NumPy code for that and @hetland has a fortran wrapped module inside octant.

I am not familiar with octant's iso functionality. That is why I am raising the questions:

1) Is it worth creating a standalone module out of iso.f? 2) Would it be possible to adapt iso.f to ugrid? 3) How difficult would it be to translate that to Cython or Numba?

hetland commented 8 years ago

1) Is it worth creating a standalone module out of iso.f?

Perhaps. I have a number of ways to do this in octant, and have never figured out the best way.

2) Would it be possible to adapt iso.f to ugrid?

Absolutely. iso.f does not care at all about the horizontal coordinate. Well, actually, it can be used to care about that, in that it creates an iso surface as a projection in the leftmost axis. In a structured grid, it is possible to do this in x, or y, but doing it in only z is the common use case, and is directly adaptable for ugrids.

3) How difficult would it be to translate that to Cython or Numba?

I expect this would be straightforward. The fortran code in iso.f is very simple (at least for linear interpolation, splines are slightly more complex). I think that a cython routine would be the right way to go here, actually. But first we need to make some design decisions:

1) linear interpolation, splines or a choice based on a kwarg. I would prefer a choice, with linear interpolation as default.

2) what to do when the field is monotonic. I would suggest the shallowest (i.e., highest index) isosurface is the right one to return. If the bottom is desired, you can always flip the vertical axis upside down. We should not require monotonicity, but we then need to be clear and consistent about what is actually returned.

Note that isosurfaces should be the preferred way to make all vertical slices. I.e., an isoslice of z=-50m should be the way to get properties at a specific depth. No need to have two different codes for slices in z, density, etc.

On Sat, Sep 19, 2015 at 1:44 PM, Filipe notifications@github.com wrote:

Isosurfaces is one of the methods raised here https://github.com/ioos/APIRUS/issues/1#issuecomment-141058183.

I have a very inefficient NumPy code for that and @hetland https://github.com/hetland has a fortran wrapped module https://github.com/hetland/octant/blob/68f51005480c8cb9f83098b0f751b4e8086c11b5/octant/src/iso.f inside octant.

I am not familiar with octant's iso functionality. That is why I am raising the questions:

1) Is it worth creating a standalone module out of iso.f? 2) Would it be possible to adapt iso.f to ugrid? 3) How difficult would it be to translate that to Cython or Numba?

— Reply to this email directly or view it on GitHub https://github.com/ioos/APIRUS/issues/6.

Prof. Rob Hetland Texas A&M Univ. – Dept. of Oceanography http://pong.tamu.edu/~rob

ocefpaf commented 8 years ago

1) linear interpolation, splines or a choice based on a kwarg. I would prefer a choice, with linear interpolation as default.

Agreed.

2) what to do when the field is monotonic. I would suggest the shallowest (i.e., highest index) isosurface is the right one to return. If the bottom is desired, you can always flip the vertical axis upside down. We should not require monotonicity, but we then need to be clear and consistent about what is actually returned.

No opinion here. I guess I never used this enough to develop one.

3) Note that isosurfaces should be the preferred way to make all vertical slices. I.e., an isoslice of z=-50m should be the way to get properties at a specific depth. No need to have two different codes for slices in z, density, etc.

Yes! That is the main reason I opened this issue and asked those questions. Based on your answer I think we are OK to lay a plan out:

I can help with the first two.

hetland commented 8 years ago

This is a good plan. I can help with 3.. as soon as I learn Cython...

On Fri, Sep 25, 2015 at 10:34 AM, Filipe notifications@github.com wrote:

1) linear interpolation, splines or a choice based on a kwarg. I would prefer a choice, with linear interpolation as default.

Agreed.

2) what to do when the field is monotonic. I would suggest the shallowest (i.e., highest index) isosurface is the right one to return. If the bottom is desired, you can always flip the vertical axis upside down. We should not require monotonicity, but we then need to be clear and consistent about what is actually returned.

No opinion here. I guess I never used this enough to develop one.

3) Note that isosurfaces should be the preferred way to make all vertical slices. I.e., an isoslice of z=-50m should be the way to get properties at a specific depth. No need to have two different codes for slices in z, density, etc.

Yes! That is the main reason I opened this issue and asked those questions. Based on your answer I think we are OK to lay a plan out:

  • Create a standalone module, Fortran wrapped, for isosurface.
  • Create a few use cases with all the grids/variables we want to see sliced.
  • Cythonize it.

I can help with the first two.

— Reply to this email directly or view it on GitHub https://github.com/ioos/APIRUS/issues/6#issuecomment-143253667.

Prof. Rob Hetland Texas A&M Univ. – Dept. of Oceanography http://pong.tamu.edu/~rob

ChrisBarker-NOAA commented 8 years ago

"""

Does "Cythonize it" mean replace the fortran with Cython, or call the fortran with Cython?

I'd rather not see a fortran dependency in this if we can avoid it.... still kindo f a pain to get the compiler all set up and working on all systems (particularly Windows!)

@hetland: I can help with the Cython. Great tool to learn.

Maybe you'll fidn this useful:

http://uwpce-pythoncert.github.io/SystemDevelopment2015/extensions.html#extensions

particularly:

https://github.com/UWPCE-PythonCert/SystemDevelopment2015/tree/master/Examples/week-08-extensions/cython/AGC_example

(note that in that example, fortran is faster than C or Cython -- but I think I know how to get C as fast - just haven't taken the time to try it out)

ocefpaf commented 8 years ago

Does "Cythonize it" mean replace the fortran with Cython, or call the fortran with Cython?

Replace the Fortran with Cython.

pain to get the compiler all set up and working on all systems (particularly Windows!)

Not with conda and mingwpy. We already have octant, with the iso.f extension, working on Windows-32. I still need to get my hands on Windows-64 to fix a failure there.

ChrisBarker-NOAA commented 8 years ago

""" Not with conda and mingwpy. We already have octant, with the iso.f extension, working on Windows-32. I still need to get my hands on Windows-64 to fix a failure there. """

I thought Win64 was still a nightmare for Fortran.... But anyway, I'm all for the Cytoninzing!

ocefpaf commented 8 years ago

The issues and TODOs are in https://github.com/pyoceans/iso.

ocefpaf commented 8 years ago

See it in action: http://nbviewer.ipython.org/gist/ocefpaf/6e9a61c8dcaa7e1523ee

(Closing this issue, but not the discuss. Let's talk at: https://github.com/pyoceans/iso.)

rsignell-usgs commented 8 years ago

See it in action: http://nbviewer.ipython.org/gist/ocefpaf/6e9a61c8dcaa7e1523ee

Super cool!!!!

ChrisBarker-NOAA commented 8 years ago

very cool. Damn you're productive! :-)

ocefpaf commented 8 years ago

I can't take credit for that one I just broke @hetland's octant apart and put it together leaving some pieces out.

(My favorite hobby, but I usually do that with electronics :wink: )