brmather / pycurious Goto Github PK

I am running a windowed data (200 by 200) km with about 40,000 points on the centroid code and I got this error:

ValueError Traceback (most recent call last)
in
8
9 window_size = 200e3
---> 10 subgrid = grid.subgrid(window_size, xpt, ypt)

~\Anaconda3\lib\site-packages\pycurious\grid.py in subgrid(self, window, xc, yc)
149 raise ValueError(
150 "Window size {} at centroid {} exceeds the data range".format(
--> 151 window, (xc, yc)
152 )
153 )

ValueError: Window size 200000.0 at centroid (600500.0, 100500.0) exceeds the data range

Part of openjournals/joss-reviews#1544.

Methods CurieGrid.radial_spectrum, CurieGrid.radial_spectrum_log and CurieGrid.azimuthal_spectrum perform almost the same task. They differ only on a few lines. It's a good practice to avoid repeating code. It makes maintenance and testing a lot easier.

Would you consider creating a private method for CurieGrid that performs common tasks for these three methods, and refactor them to just call this new private method?

Part of openjournals/joss-reviews#1544.

I'm a little bit confused about the scale argument of CurieGrid.radial_spectrum:

Is it only used for units conversion? If that's the case, it's not a good choice to let users to change unit conversions constants.

By the way, is there a special reason why converting units from meters (grid coordinates) to km (wavenumbers and estimated depths)?

I was not able to find the 'grid.optimize_routine' function w.r.t the Tanaka part of the code. The objective is to be able to view a curie depth plot for my entire study region using Tanaka's method. Can you help, Dr. Mather?

Part of openjournals/joss-reviews#1544.

Because pycurious relies on processing the magnetic data in order to compute the CDP:

• read/load magnetic grid,
• select a subgrid,
• compute FFT and its radial average,
• apply tappering and
• estimate depths to magnetic sources,

I consider relevant to extend tests for the core computations, just to make sure they work as expected and to prevent breaking anything on future changes.
The current test functions only test if the methods don't fail when run. That's important, but I think is also mandatory to check if they are working as expected.
I don't intent you to achieve a 100% coverage, but just add tests that assert if the result is close to the expected one.

Besides individual tests for each method, would be nice to add a synthetic magnetic data (with known depths of the magnetic source) and check if Tanaka and Bouligand methods recover the original depths.

What do you think? Do you agree?

I am having a strange fractal exponent (8.34) when I apply emag2_v3 to bouligand approach. The window size is 200km by 200km.
Please help out with suggestions.
Cheers

hello, I have a question about your code to compute radially averaged power spectrum.
In grid/curie.py, you compute power spectrum in function radial_spectrum

FT = np.abs(np.fft.fft2(data*vtaper))

But as I known, it shoud be FT = np.abs(np.fft.fft2(data*vtaper))**2, the square of absolute of fft, isn't it?

Nice package you have here! I didn't know it. I wrote some messy code for Curie Depth Point calculation when I was undergraduate, but it's old dated, not well written nor documented (my first baby steps on Python).
It implements the Tanaka method and a new I've been experimenting with: fitting the hole exponential curve with a Levenberg-Marquardt algorithm.
Besides I wrote an ugly GUI to manually select windows and seeing radial averaged power spectrum live.
I think I haven't uploaded the final version, but if you are interested I can do it.

Nevermind, your package looks really good!
I'll recommend yours if someone ask me to compute CDP!

Hello!

My name is Reed, and I'm a software engineering researcher at Sandia National Laboratories in the US. I've created an issue on your repo just to ask a quick question. If you don't have time or don't care to respond, feel free to ignore me and/or delete this issue.

Where I work, we have a very diverse ecosystem of cutting-edge research codes spanning every discipline you could imagine. I'm part of our software engineering research department, and it's our job to keep that ecosystem robust and healthy. Part of that means helping our scientists to adopt good software practices. Right now, my mind has been on software versioning/release schemes (e.g. semantic versioning).

In order to build a case for/against getting my people on-board with the practice, I figured I should ask people who already use versioning to release their software to see what they think. So I gathered up a list of scientific software repositories on GitHub, then I selected those that tracked versioned releases, that were reasonably active, etc. Finally, I picked a handful of those repos and decided to reach out to them. You were on that list.

Anyway, here's the question:

What do you believe are the benefits (or drawbacks) of having versioned releases of your software (i.e. 1.0.0, 1.1.0, 2.0.0...)? When should someone start thinking about versioning/releasing their code?

Just a sentence or two, that's all I need. For context, imagine the preceding sentence is this: "But don't just take my word for it, just listen to what this accomplished computational geophysicist has to say!"

Thank you so much!

Reed Milewicz
[email protected]

Part of openjournals/joss-reviews#1544.

There are two unused variables on tanaka1999():

Are they not necessary or are missing on any following line?

Part of openjournals/joss-reviews#1544.

When checking if the node spacing is equal on both northing and easting directions, a maximum difference of 1 meter has been hardcoded:

Would be nice to use np.allclose(dx, dy) instead of hardcoding a minumum difference that could raise some unwanted results. For example, if you have a very dense survey with one measurement every 10m and dx is around 1m different from dy, you'll finally get a significant error on the wavenumbers and therefore on the estimated depths even that error won't be reflected on the statistical errors for these depths.

Do you agree?

Besides, wouldn't be better to raise an error instead of a warning? Using non equal spaced grids is acceptable for further computations?

On Ex5-Mapping-Curie-depth-EMAG2.ipynb an error is raised when running the following cell:

grid = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)

window_size = 400e3

# centroid spacing of 50 km
xc_list, yc_list = grid.create_centroid_list(window_size, spacingX=50e3, spacingY=50e3)
print("number of centroids = {}".format(len(xc_list)))


beta, zt, dz, C = grid.optimise_routine(window_size, xc_list, yc_list)


>>> ValueError: node spacing should be identical (1000.5265929436546, 1000.5885815185403)

Seems that when creating the grid, the nodes are not equally spaced between two axis.

Part of openjournals/joss-reviews#1544.

Sometimes the docstrings refer to Zor as the bottom of the magnetic source. Nevertheless, when computing Zb, ComputeTanaka() applies the following:

Is the docstring mistaken or I'm missing something here?

References

• Tanaka A., Okubo Y., Matsubayashi O. (1999), Curie point depth based on spectrum analysis of the magnetic anomaly data in East and Southeast Asia, Tectonophysics, doi: 10.1016/S0040-1951(99)00072-4

Part of openjournals/joss-reviews#1544.

On CurieGrid docstring, would be nice to specify what x and y stands for.
Some users may be used to packages that uses y for easting and x for northing, while others may refer to a classical Cartesian axes, where x is usually assigned to the horizontal (easting) coordinate and y to the vertical one (northing).

Also, would be nice to specify the expected shape of the grid 2d-array: (ny, nx).

Part of openjournals/joss-reviews#1544.

CuriedGrid has the following methods: radial_spectrum and radial_spectrum_log.
Their names suggest the second one computes the radial spectrum and applies the np.log, while the first one does not. Nevertheless, the difference between both of them does not rely on the np.log but on the np.sqrt:

CurieGrid.radial_spectrum

CurieGrid.radial_spectrum_log

Moreover, I don't quite understand what's the goal of CurieGrid.azimuthal_spectrum. That could be related to some lack of knowledge on the subject from my part.

I think extending their docstring would help to a better understanding of what they are doing. Maybe including some math and references would be nice. Renaming the methods in order to make very clear what they do it's also a good approach.

Am I understanding something wrong here or do you agree with my point of view?

Part of openjournals/joss-reviews#1544.

When extracting a subgrid, an error is raised if its center falls outside the data range:

That's a nice thing to check! I would add another check: if the window is too big it may reach regions beyond the data range, what may cause errors afterwards but would be hard to track where they were originated.
What do you think?

Part of openjournals/joss-reviews#1544.

When trying to install pycurious with pip on Python 3 (Anaconda) I get the following error:

$ pip install pycurious
Collecting pycurious
  Using cached https://files.pythonhosted.org/packages/2a/f4/063241be23baf131ef7d6967d26d353fa9a00eb4fe0751db9b991df944e4/pycurious-0.4.tar.gz
Requirement already satisfied: numpy in ./.anaconda3/envs/pycurious/lib/python3.7/site-packages (from pycurious) (1.16.4)
Requirement already satisfied: scipy>=0.15.0 in ./.anaconda3/envs/pycurious/lib/python3.7/site-packages (from pycurious) (1.3.0)
Requirement already satisfied: Cython>=0.25.2 in ./.anaconda3/envs/pycurious/lib/python3.7/site-packages (from pycurious) (0.29.12)
Building wheels for collected packages: pycurious
  Building wheel for pycurious (setup.py) ... error
  ERROR: Complete output from command /home/santi/.anaconda3/envs/pycurious/bin/python -u -c 'import setuptools, tokenize;__file__='"'"'/tmp/pip-install-1e1fzccz/pycurious/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' bdist_wheel -d /tmp/pip-wheel-dyyg25o2 --python-tag cp37:
  ERROR: running bdist_wheel
  running build
  running build_py
  creating build
  creating build/lib.linux-x86_64-3.7
  creating build/lib.linux-x86_64-3.7/pycurious
  copying pycurious/download.py -> build/lib.linux-x86_64-3.7/pycurious
  copying pycurious/grid.py -> build/lib.linux-x86_64-3.7/pycurious
  copying pycurious/__init__.py -> build/lib.linux-x86_64-3.7/pycurious
  copying pycurious/documentation.py -> build/lib.linux-x86_64-3.7/pycurious
  copying pycurious/mapping.py -> build/lib.linux-x86_64-3.7/pycurious
  copying pycurious/optimise.py -> build/lib.linux-x86_64-3.7/pycurious
  creating build/lib.linux-x86_64-3.7/pycurious/Examples
  creating build/lib.linux-x86_64-3.7/pycurious/Examples/data
  copying pycurious/Examples/data/test_mag_data.txt -> build/lib.linux-x86_64-3.7/pycurious/Examples/data
  copying pycurious/Examples/0-StartHere.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples
  creating build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks
  creating build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Tanaka
  copying pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Tanaka
  copying pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Tanaka
  copying pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Tanaka
  creating build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Bouligand
  copying pycurious/Examples/Notebooks/Bouligand/Ex4-Spatial-variation-of-Curie-depth.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Bouligand
  copying pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Bouligand
  copying pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Bouligand
  copying pycurious/Examples/Notebooks/Bouligand/Ex5-Mapping-Curie-depth-EMAG2.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Bouligand
  copying pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb -> build/lib.linux-x86_64-3.7/pycurious/Examples/Notebooks/Bouligand
  running build_ext
  building 'pycurious.radon' extension
  creating build/temp.linux-x86_64-3.7
  creating build/temp.linux-x86_64-3.7/src
  gcc -pthread -B /home/santi/.anaconda3/envs/pycurious/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/santi/.anaconda3/envs/pycurious/lib/python3.7/site-packages/numpy/core/include -Isrc -I/home/santi/.anaconda3/envs/pycurious/include/python3.7m -c src/radon.c -o build/temp.linux-x86_64-3.7/src/radon.o -std=c99
  In file included from /home/santi/.anaconda3/envs/pycurious/lib/python3.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1822,
                   from /home/santi/.anaconda3/envs/pycurious/lib/python3.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,
                   from /home/santi/.anaconda3/envs/pycurious/lib/python3.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                   from src/radon.c:625:
  /home/santi/.anaconda3/envs/pycurious/lib/python3.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning "Using deprecated NumPy API, disable it with " "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
     17 | #warning "Using deprecated NumPy API, disable it with " \
        |  ^~~~~~~
  gcc -pthread -B /home/santi/.anaconda3/envs/pycurious/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/santi/.anaconda3/envs/pycurious/lib/python3.7/site-packages/numpy/core/include -Isrc -I/home/santi/.anaconda3/envs/pycurious/include/python3.7m -c src/cradon.c -o build/temp.linux-x86_64-3.7/src/cradon.o -std=c99
  src/cradon.c: In function ‘radon2d’:
  src/cradon.c:228:9: warning: returning ‘void *’ from a function with return type ‘int’ makes integer from pointer without a cast [-Wint-conversion]
    228 |         import_array();  // to use PyArray_SimpleNewFromData
        |         ^~~~~~~~~~~~
  gcc -pthread -shared -B /home/santi/.anaconda3/envs/pycurious/compiler_compat -L/home/santi/.anaconda3/envs/pycurious/lib -Wl,-rpath=/home/santi/.anaconda3/envs/pycurious/lib -Wl,--no-as-needed -Wl,--sysroot=/ build/temp.linux-x86_64-3.7/src/radon.o build/temp.linux-x86_64-3.7/src/cradon.o -Lpycurious -o build/lib.linux-x86_64-3.7/pycurious/radon.cpython-37m-x86_64-linux-gnu.so
  /home/santi/.anaconda3/envs/pycurious/compiler_compat/ld: build/temp.linux-x86_64-3.7/src/radon.o: unable to initialize decompress status for section .debug_info
  /home/santi/.anaconda3/envs/pycurious/compiler_compat/ld: build/temp.linux-x86_64-3.7/src/radon.o: unable to initialize decompress status for section .debug_info
  /home/santi/.anaconda3/envs/pycurious/compiler_compat/ld: build/temp.linux-x86_64-3.7/src/radon.o: unable to initialize decompress status for section .debug_info
  /home/santi/.anaconda3/envs/pycurious/compiler_compat/ld: build/temp.linux-x86_64-3.7/src/radon.o: unable to initialize decompress status for section .debug_info
  build/temp.linux-x86_64-3.7/src/radon.o: file not recognized: file format not recognized
  collect2: error: ld returned 1 exit status
  error: command 'gcc' failed with exit status 1
  ----------------------------------------
  ERROR: Failed building wheel for pycurious
  Running setup.py clean for pycurious
Failed to build pycurious

I don't think the two raised warnings are causing the error. I'm more worried about the file format not recognized error when reading radon.o.

Do you know what am I doing wrong or how to solve this?

When solving [Ex5-Mapping-Curie-depth-EMAG2.ipynb] through Binder, at the interpolation onto grid projection step, I got this error:

ModuleNotFoundError Traceback (most recent call last)
in
7
8 # extent on the sphere (WGS84)
----> 9 extent_sphere = mapping.convert_extent(extent_grid, epsg_in=2157, epsg_out=4326)
10
11 # map extent - should be narrower than extent_sphere

/usr/lib/python3.6/site-packages/pycurious/mapping.py in convert_extent(extent_in, epsg_in, epsg_out)
103 xi = [xmin, xmin, xmax, xmax]
104 yi = [ymin, ymax, ymin, ymax]
--> 105 xo, yo = transform_coordinates(xi, yi, epsg_in, epsg_out)
106 extent_out = [min(xo), max(xo), min(yo), max(yo)]
107 return extent_out

/usr/lib/python3.6/site-packages/pycurious/mapping.py in transform_coordinates(x, y, epsg_in, epsg_out)
77 y coordinates projected in epsg_out
78 """
---> 79 import pyproj
80
81 proj_in = pyproj.Proj("+init=EPSG:" + str(epsg_in))

ModuleNotFoundError: No module named 'pyproj'

Could you please debug this?

Part of openjournals/joss-reviews#1544.

When computing the error for the depth to the bottom of the magnetic source, the errors of Zor and Ztr should be added instead of subtracted:

 dZb = 2.0*dZor + dZtr

Part of openjournals/joss-reviews#1544.

Is there a special reason why data is copied on the following line?

Remember that the copy() method allocates memory, what would increase the memory consumption if data is too big and require more computation time.