Closed gerritholl closed 4 years ago
gerritholl
commented
4 years ago Some more information on why I suspect resampling is to blame
ir_123.tif) the problem is identicalgerritholl
commented
4 years ago I welcome ideas on how I can find out more clues.
djhoese
commented
4 years ago Thoughts:
PYTROLL_CHUNK_SIZE=5568).gerritholl
commented
4 years ago When I divide each component of area_extent by 10, there is no more segmentation fault, instead there is an AssertionError in Scene._slide_data:
Traceback (most recent call last):
File "mwe41.py", line 17, in <module>
nsc = sc.resample(area)
File "/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/satpy/scene.py", line 1111, in resample
self._resampled_scene(new_scn, destination, resampler=resampler,
File "/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/satpy/scene.py", line 1055, in _resampled_scene
dataset = self._slice_data(source_area, (slice_x, slice_y), dataset)
File "/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/satpy/scene.py", line 993, in _slice_data
assert ('y', source_area.y_size) in dataset.sizes.items()
AssertionErrorThe full log for this situation:
Code leading to the error log:
from satpy import Scene
from satpy.utils import debug_on
from pyresample import create_area_def
debug_on()
from glob import glob
area = create_area_def("fci_0deg_2km",
{"proj": "geos", "lon_0": 0, "h": 35786400, "x_0": 0, "y_0": 0, "ellps": "WGS84"},
units="m",
shape=(5568, 5568),
area_extent=(-5567999.994206558/10, -5567999.994200589/10,
5567999.994206558/10, 5567999.994200589/10))
sc = Scene(glob("/media/nas/x21308/2020_04_MTG_unofficial_Testdata/20130804_RC72/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--L2P-NC4E_C_EU*.nc"), reader="fci_l1c_fdhsi")
print("doing load + compute")
sc.load(["ir_123"])
sc["ir_123"].compute()
print("doing resample + compute")
nsc = sc.resample(area)
nsc["ir_123"].compute()
print("saving")
nsc.save_dataset("ir_123", "/tmp/ir_123.tif")
pnuu
commented
4 years ago Did you also divide the shape (width and height) by 10?
gerritholl
commented
4 years ago Results are identical with Python 3.7.
I get the same AssertionError if I also divide the shape (both width and height) by 10 (I initially didn't).
I wonder if this message (the last message before the AssertionError in the detailed log, and also included in the log preceding the segmentation fault) provides a clue:
[DEBUG: 2020-05-06 16:06:36 : pyresample.geometry] Projections for data and slice areas are identical: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1],ID["EPSG",7030]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Geostationary Satellite (Sweep Y)"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Satellite Height",35786400,LENGTHUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]Going into the debugger at the AssertionError, the shape of dataset is apparently (0, 558), with the value of dataset.sizes that is being asserted:
ipdb> p dataset.sizes
Frozen({'y': 0, 'x': 558})The source_area has a negative number of rows:
ipdb> p source_area
Area ID: some_area_name
Description: On-the-fly area
Projection ID: geosfci
Projection: {'ellps': 'WGS84', 'h': '35786400', 'lon_0': '0', 'no_defs': 'None', 'proj': 'geos', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 558
Number of rows: -556
Area extent: (-557999.9994, -555999.9994, 557999.9994, 555999.9994)Ah, that's a known issue: https://github.com/pytroll/satpy/issues/930
pnuu
commented
4 years ago And one workaround is to use nsc = sc.resample(area, reduce_data=False).
Maybe I should dedicate the rest of the PCW to complete the PR I started in October.
pnuu
commented
4 years ago I'm continuing my old work in https://github.com/pytroll/satpy/pull/932 which should fix this.
gerritholl
commented
4 years ago I am not sure if the AssertionError and the segmentation fault have the same cause. #275 solves the former but not the latter. See my comments there for more information.
pnuu
commented
4 years ago I started to investigate this.
The area definition @gerritholl has defined has these extents
(-5567999.994206558, -5567999.994200589, 5567999.994206558, 5567999.994200589))
while the extents in the data are
(-5567999.994206558, 5567999.994206558, 5567999.994200589, -5567999.994200589)
and when the area def is flipped to the same orientation Gerrit has them (North up) with attrs_area[::-1] they become
(-5567999.994206558, -5569999.994198508, 5567999.994200588, 5565999.9942086395)
which are a bit different from the extents Gerrit has defined. The second and fourth values differ by almost one resolution unit:
area.area_extent[1] - area.resolution[1] - attrs_area_flipped.area_extent[1] -> 1.862645149230957e-09
and
area.resolution[1] - attrs_area_flipped.resolution[1] -> -1.071839506039396e-09
The area definition @gerritholl made seems to be valid, apart from being slightly different from the area def the data has, and the dimensions match:
ipdb> (area.area_extent[2] - area.area_extent[0]) / area.resolution[0]
5568.0
ipdb> (area.area_extent[3] - area.area_extent[1]) / area.resolution[1]
5568.0The flipped area from the attributes works when resampling:
attrs_area = glbl['ir_123'].attrs['area']
# flip the area so that North is up
attrs_area_flipped = attrs_area[::-1, :]
lcl = glbl.resample(attrs_area_flipped)
No idea why the area @gerritholl defined doesn't work... :thinking: :man_shrugging:
pnuu
commented
4 years ago Oh! In the area definition of the data, both resolutions are equal:
attrs_area_flipped.resolution[0] - attrs_area_flipped.resolution[1] -> 0.0
whereas in @gerritholl's there is a small difference:
area.resolution[0] - area.resolution[1] -> 2.143906385754235e-09
which indicates that the extents do not cover the same distance in both x and y directions. And yes, there is a difference of 11.9 micrometers:
(area.area_extent[2] - area.area_extent[0]) - (area.area_extent[3] - area. area_extent[1]) -> 1.1937692761421204e-05.
Not sure why, but for some reason pykdtree doesn't like this. @mraspaud any ideas?
pnuu
commented
4 years ago And this is why @gerritholl and I think it is pykdtree, the output in dmesg after a segfault:
[298457.507131] python[13687]: segfault at 560433c31380 ip 00007fa2a03a6653 sp 00007fa29378d068 error 4 in kdtree.cpython-37m-x86_64-linux-gnu.so[7fa2a0397000+11000]
pnuu
commented
4 years ago Some more debugging. The area definition @gerritholl defined ends up being with a shape (1, 5568) in pyresample.kd_tree.XArrayResamplerNN. The flipped original has a full shape of (5568, 5568).
gerritholl
commented
4 years ago It turns out that despite merging #277 this problem is not solved. Although I don't get a segfault storing ir_123 anymore, I do get one storing day_microphysics:
from satpy import Scene
from glob import glob
from satpy.utils import debug_on
debug_on()
path_to_testdata = "/media/nas/x21308/2020_04_MTG_unofficial_Testdata/20130804_RC72/"
scn = Scene(filenames=glob(path_to_testdata + "/*BODY*.nc"),
reader=['fci_l1c_fdhsi'])
scn.load(['day_microphysics'])
sc = scn.resample('fci_0deg_2km', reduce_data=True)
sc["day_microphysics"].persist()
sc.save_dataset('day_microphysics', filename=f"/tmp/test.tiff")Running with python -X dev gives at the very end:
The full log is too long to fit in a comment, but available here: log.txt
Andrea reports that he still gets a segmentation fault even storing ir_123. All versions are latest master:
pyresample 1.15.0+192.g1c2b794 pypi_0 pypi
pyspectral 0.9.6.dev17+g7d30889 pypi_0 pypi
satpy 0.21.1.dev290+ga761d23d pypi_0 pypi
trollimage 1.12.0+15.g1f63856 pypi_0 pypiAnd indeed dmesg has
[11861417.179337] python[31849]: segfault at 558f4a6c7e20 ip 00007f3758471033 sp 00007f37337fd098 error 4 in kdtree.cpython-38-x86_64-linux-gnu.so[7f3758460000+13000]Segfault confirmed, even with ir_123 channel data. I'm having a look.
pnuu
commented
4 years ago The fix I made in #277 are not in Pyresample master. Re-adding those changes locally makes the segfault go away for 'ir_123', but does crash for glbl.load(['ir_123', 'natural_color', 'airmass']). Inspecting further.
pnuu
commented
4 years ago Out of these three, it's the 'natural_color' composite that triggers the segfault in pykdtree, other two work just fine.
pnuu
commented
4 years ago Oh, 'vis_08_pixel_quality' is being resampled, I bet this causes the crash.
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Computing kd-tree parameters
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Resampling pixel_quality
[WARNING: 2020-06-04 08:14:07 : pyresample.kd_tree] Fill value incompatible with integer data using 255 instead.
/home/lahtinep/Software/pytroll/pytroll_packages/pyresample/pyresample/kd_tree.py:1235: PerformanceWarning: Increasing number of chunks by factor of 11
dtype=new_data.dtype, concatenate=True)
[DEBUG: 2020-06-04 08:14:07 : satpy.scene] Resampling DatasetID(name='vis_08', wavelength=(0.815, 0.865, 0.915), resolution=1000, polarization=None, calibration='reflectance', level=None, modifiers=())
/home/lahtinep/Software/miniconda3/envs/pytroll/lib/python3.7/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj_string = self.to_proj4()
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Check if ./resample_lut-8f7af1b20a2425c92db67dc8fdf8d46ca9c02df6.npz exists
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Computing kd-tree parameters
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Resampling None
/home/lahtinep/Software/pytroll/pytroll_packages/pyresample/pyresample/kd_tree.py:1235: PerformanceWarning: Increasing number of chunks by factor of 11
dtype=new_data.dtype, concatenate=True)
[DEBUG: 2020-06-04 08:14:07 : satpy.scene] Resampling DatasetID(name='vis_08_pixel_quality', wavelength=None, resolution=1000, polarization=None, calibration=None, level=None, modifiers=())
/home/lahtinep/Software/miniconda3/envs/pytroll/lib/python3.7/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj_string = self.to_proj4()
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Check if ./resample_lut-8f7af1b20a2425c92db67dc8fdf8d46ca9c02df6.npz exists
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Computing kd-tree parameters
[DEBUG: 2020-06-04 08:14:07 : satpy.resample] Resampling pixel_quality
[WARNING: 2020-06-04 08:14:07 : pyresample.kd_tree] Fill value incompatible with integer data using 255 instead.
/home/lahtinep/Software/pytroll/pytroll_packages/pyresample/pyresample/kd_tree.py:1235: PerformanceWarning: Increasing number of chunks by factor of 11
dtype=new_data.dtype, concatenate=True)And when trying to load that dataset:
glbl.load(['vis_08_pixel_quality'])
*** KeyError: 'data/vis_08_pixel_quality/measured/effective_radiance/shape'It is shown in the output of glbl.available_dataset_names(). Huh?
pnuu
commented
4 years ago Loading, resampling and saving 'vis_08' triggers the segfault, nothing else is needed.
pnuu
commented
4 years ago Trying to resample any of the vis_* or nir_* channels triggers the segfault. All wv_* and ir_* channels are fine.
gerritholl
commented
4 years ago It seems that some of the magic I'm doing in https://github.com/pytroll/satpy/pull/1177 is causing problems. Each channel data variable in the FCI data is located in its own group. Each of these groups has a data variable pixel_quality, which is referred to in the ancillary_variables data variable attribute for the corresponding channel. Because we can't load 16 data variables that are all called pixel_quality and we don't have another way to properly disambugate them before https://github.com/pytroll/satpy/pull/1088 (dynamic dataset IDs) is merged, the FCI reader rewrites the contents of the pixel_quality ancillary variable when channel data are being read, to prepend the name of the channel. But these data variables are not actually in the data under that name, and can't be loaded explicitly: they can only be taken in by reading the corresponding channel.
This is not a pretty solution but it's a temporary one we settled upon until we can use dynamic dataset IDs to define a field for the pixel_quality data variables indicating the parent datasetid this pixel quality dataset refers to.
I wonder if this magic somehow interferes with the resampling, or with the implicit loading of channels as composite dependencies? But that wouldn't explain why it affects some channels but not others, since this magic should equally affect all channels.
gerritholl
commented
4 years ago Note: vis and nir have a 1 km resolution, wv and ir have a 2 km resolution. The equivalent test of "resampling" each to the own area would be "resampling" vis and nir channels to fci_0deg_1km, and wv and ir channels to fci_0deg_2km.
pnuu
commented
4 years ago Also resampling to 1 km grid segfaults. And the nearest neighbour resampling should in any case work, no matter the resolution.
gerritholl
commented
4 years ago It seems that when using https://github.com/pytroll/satpy/pull/1230 the segmentation fault does not occur (I have no clue why).
pnuu
commented
4 years ago Also #285 fixes the segmentation fault.
Code Sample, a minimal, complete, and verifiable piece of code
Problem description
The code snippet results in a segmentation fault (see details below).
Expected Output
I expect an image with resampled to the indicated area. If there is a problem with my area definition, I expect an informative error message and a graceful exit.
Actual Result, Traceback if applicable
Output of python -X dev mwe41.py
```Versions of Python, package at hand and relevant dependencies