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

Add baseline for GeoColor composite including FCI, AHI and ABI recipes #2557

Merged
merged 54 commits into from
Dec 15, 2023

Conversation

strandgren
Copy link
Collaborator

@strandgren strandgren commented Aug 24, 2023

This PR adds the baseline for the GeoColor composite blend developed by Miller et al. 2020. The GeoColor implementation in this PR follows the general approach described in Miller et al., but uses equivalent or similar Satpy functionalities, methods and composites where available. Some notable differences are:

  • The green correction used for the true color layer is based on the hybrid green method used in Miller et al., but including the pixel-level scaling by NDVI available in Satpy.
  • The true color layer for ABI is currently using the already available true_color composite based in the green channel simulation in Satpy, with no use of the green channel simulation approach and hybrid green correction used in Miller et al. (2020).
  • We do not account for cloud height information in the IR channel when doing the Rayleigh correction. Currently, there is also no blend with an uncorrected image to improve the imagery at high sunz angles (does not appear to be needed though).
  • The night-time layer used here is simply the NASA Black Marble geotiff image. In Miller et al. the use an elevation-enhanced nightscape layer with lights data from VIIRS. They also use an orange color for the city light to better replicate the appearance sodium lamps commonly used in the U.S.

The PR implements the missing compositors/layers needed for the GeoColor composite as well as recipes for FCI, AHI and ABI. Based on this, it should be easy to add recipes for other sensors like AMI, AGRI or even polar orbiting satellites in the future.

Below two early examples for FCI and AHI:

FCI
image

AHI
image

…rence used for the Geocolor low-level cloud layer.
… of projectables and new land-sea mask data.
@strandgren
Copy link
Collaborator Author

strandgren commented Aug 24, 2023

There are currently two aspects that need more work (apart from the tests..):

  1. The night-time low-level cloud layer requires a land-sea mask. I have created a geotiff for this using the gshhs dataset at intermediate resolution. Currently this file is stored locally, but we would need to add it to an online repository and implement automatic downloading and caching (same as for NASA Black Marble). Here I need support for adding the land-sea mask file to a suitable online repository. The file can be accessed here if anyone is interested in testing the code in this PR: https://sftp.eumetsat.int/public/file/2G3bKqSCpEW40GYNtLJ3EQ/gshhs_land_sea_mask_3km_i.tif. Sea pixel have value 0 and land pixels value 100.

  2. The spectral test for identifying low-level clouds at night is based on a single brightness temperature difference test, which can lead to many false alarms over barren surface types like Sahara and Australia (most likely because of spectral emissivity differences). Furthermore, we might have to wait until later during FCI commissioning in order to assess the quality of the low-level cloud detection for FCI and tune the detection thresholds accordingly. Nevertheless, I would be in favor of merging as is (keeping the detection thresholds from Miller et al., 2020 as default) and then work on this more in a separate PR. Ideally together with NOAA/CIRA and anyone that has the time to test the current implementation and has ideas for improvements :)

UPDATE:

  1. The global land-water mask GeoTIFF file has been uploaded to zenodo (https://zenodo.org/records/10076199) and this file is now being automatically downloaded and cached when creating a GeoColor composite.

  2. The default thresholds used for the high- and low-level cloud layer detection have been updated based on feedback from CIRA (GeoColor developers). Final thresholds for FCI will have to wait for calibrated FCI data. False alarms and missed clouds is still a problem for the low-level cloud layer. This is a flaw of the simple spectral test approach used here. An improved and more accurate low-level cloud detection is being worked on (both by CIRA and also at EUM), but will have to go into a separate PR in the future.

@djhoese djhoese added enhancement code enhancements, features, improvements component:compositors labels Aug 24, 2023
@strandgren strandgren changed the title Add baseline for GeoColor composite including FCI and AHI recipes Add baseline for GeoColor composite including FCI, AHI and ABI recipes Aug 25, 2023
@strandgren
Copy link
Collaborator Author

strandgren commented Aug 25, 2023

And here a couple of examples with ABI:

ABI Sunrise
image

ABI Sunset
image

@yukaribbba
Copy link
Contributor

I love this.

@yukaribbba
Copy link
Contributor

There's also another option for landmask -- NASA's blue marble next-gen. https://neo.gsfc.nasa.gov/archive/bluemarble/bmng/landmask/. That huge world.watermask.86400x43200.bin.gz is in about 500m resolution. The file is binary so you need an ENVI header to read and convert it to geotiff.

…add_GeoColor_composite

Conflicts:
	satpy/etc/composites/abi.yaml
	satpy/etc/composites/fci.yaml
@codecov
Copy link

codecov bot commented Sep 1, 2023

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (c5c18cf) 95.32% compared to head (8f44b31) 95.38%.
Report is 29 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #2557      +/-   ##
==========================================
+ Coverage   95.32%   95.38%   +0.06%     
==========================================
  Files         371      371              
  Lines       52441    52521      +80     
==========================================
+ Hits        49988    50097     +109     
+ Misses       2453     2424      -29     
Flag Coverage Δ
behaviourtests 4.16% <4.72%> (-0.01%) ⬇️
unittests 96.00% <100.00%> (+0.06%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@mraspaud
Copy link
Member

While I think about it, could you add a description of geocolor to this page of the wiki? https://github.com/pytroll/satpy/wiki/Built-in-composites

@strandgren
Copy link
Collaborator Author

@yukaribbba Can you please also provide the code that triggers the error?

from satpy import config, Scene, find_files_and_readers, utils
from pyresample.geometry import AreaDefinition

config.set(config_path=['D:/satpy_config'])
utils.debug_on()
files = find_files_and_readers(base_dir='C:/Users/45107/Downloads/Sat/Geo/H08_20220101_1400_FLDK', reader='ahi_hsd')

FLDK_area = AreaDefinition(area_id='fldk', description='full disk', proj_id='geos',
                           projection='+proj=geos +h=35785863 +lon_0=140.7 +a=6378137.0 +rf=298.257024882273',
                           width=5500, height=5500, area_extent=(-5499999.968, -5499999.968, 5499999.968, 5499999.968))

scn1 = Scene(filenames=files)
scn1.load(['high_cloud_alpha_2000'])
scn1 = scn1.resample(FLDK_area, resampler='nearest')
scn1.show('high_cloud_alpha_2000')

scn2 = Scene(filenames=files)
scn2.load(['low_cloud_alpha_2000'])
scn2 = scn2.resample(FLDK_area, resampler='nearest')
scn2.show('low_cloud_alpha_2000')

scn3 = Scene(filenames=files)
scn3.load(['low_cloud_alpha_2000_with_background'])
scn3 = scn3.resample(FLDK_area, resampler='nearest')
scn3.show('low_cloud_alpha_2000_with_background')

scn4 = Scene(filenames=files)
scn4.load(['night_cloud_alpha_2000_with_background'])
scn4 = scn4.resample(FLDK_area, resampler='nearest')
scn4.show('night_cloud_alpha_2000_with_background')

I resample the scenes so that they can match the same area extent with landmask/blackmarble. And It's getting worse now. Only scn1 succeeded, all other three failed.

[DEBUG: 2023-12-14 17:50:19 : satpy.readers.yaml_reader] Requested orientation for Dataset B07 is 'native' (default). No flipping is applied.
[DEBUG: 2023-12-14 17:50:19 : satpy.readers.yaml_reader] Reading ('C:\\Users\\45107\\miniforge3\\Lib\\site-packages\\satpy\\etc\\readers\\generic_image.yaml',)
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] PROJ data files are available at built-in paths.
[DEBUG: 2023-12-14 17:50:19 : satpy.readers.yaml_reader] Assigning to generic_image: [<FSFile "C:/Users/45107/Downloads/Sat/Geo/satpy_landmask.tif">]
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entering env context: <rasterio.env.Env object at 0x000001D553BE30D0>
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Starting outermost env
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] No GDAL environment exists
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] New GDAL environment <rasterio._env.GDALEnv object at 0x000001D5014CDF00> created
[DEBUG: 2023-12-14 17:50:19 : rasterio._filepath] Installing FilePath filesystem handler plugin...
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] GDAL_DATA found in environment.
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] PROJ data files are available at built-in paths.
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x000001D5014CDF00>.
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entered env context: <rasterio.env.Env object at 0x000001D553BE30D0>
[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Sharing flag: 0
[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Nodata success: 0, Nodata value: 0.000000
[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Dataset <open DatasetReader name='C:/Users/45107/Downloads/Sat/Geo/satpy_landmask.tif' mode='r'> is started.
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting env context: <rasterio.env.Env object at 0x000001D553BE30D0>
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Cleared existing <rasterio._env.GDALEnv object at 0x000001D5014CDF00> options
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Stopped GDALEnv <rasterio._env.GDALEnv object at 0x000001D5014CDF00>.
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting outermost env
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exited env context: <rasterio.env.Env object at 0x000001D553BE30D0>
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entering env context: <rasterio.env.Env object at 0x000001D502DA3910>
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Starting outermost env
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] No GDAL environment exists
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] New GDAL environment <rasterio._env.GDALEnv object at 0x000001D502DA8520> created
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] GDAL_DATA found in environment.
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] PROJ data files are available at built-in paths.
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x000001D502DA8520>.
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entered env context: <rasterio.env.Env object at 0x000001D502DA3910>
[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Sharing flag: 0
[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Nodata success: 0, Nodata value: 0.000000
[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Dataset <open DatasetReader name='C:/Users/45107/Downloads/Sat/Geo/satpy_landmask.tif' mode='r'> is started.
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting env context: <rasterio.env.Env object at 0x000001D502DA3910>
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Cleared existing <rasterio._env.GDALEnv object at 0x000001D502DA8520> options
[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Stopped GDALEnv <rasterio._env.GDALEnv object at 0x000001D502DA8520>.
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting outermost env
[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exited env context: <rasterio.env.Env object at 0x000001D502DA3910>
C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\dataset.py:278: UserWarning: The specified chunks separate the stored chunks along dimension "x" starting at index 4096. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
[DEBUG: 2023-12-14 17:50:19 : satpy.composites.config_loader] Looking for composites config file images.yaml
[DEBUG: 2023-12-14 17:50:19 : satpy.composites.config_loader] No composite config found called images.yaml
[DEBUG: 2023-12-14 17:50:19 : satpy.readers.generic_image] Reading 'image.'
[WARNING: 2023-12-14 17:50:19 : satpy.readers.generic_image] cannot convert float NaN to integer
[DEBUG: 2023-12-14 17:50:19 : satpy.composites] Not all areas are the same in 'low_cloud_alpha_2000'
[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Delaying generation of DataID(name='low_cloud_alpha_2000') because of incompatible areas
[WARNING: 2023-12-14 17:50:19 : satpy.scene] The following datasets were not created and may require resampling to be generated: DataID(name='low_cloud_alpha_2000')
[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Unloading dataset: DataID(name='B07', wavelength=WavelengthRange(min=3.7, central=3.9, max=4.1, unit='µm'), resolution=2000, calibration=<2>, modifiers=())
[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Resampling DataID(name='B13', wavelength=WavelengthRange(min=10.2, central=10.4, max=10.6, unit='µm'), resolution=2000, calibration=<2>, modifiers=())
[DEBUG: 2023-12-14 17:50:19 : pyresample.geometry] Projections for data and slice areas are identical: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378137,298.257024882273,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Geostationary Satellite (Sweep Y)"],PARAMETER["Longitude of natural origin",140.7,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Satellite Height",35785863,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]]]]
[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Resampling DataID(name='_low_cloud_alpha_2000_dep_0', resolution=2000)
[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Resampling DataID(name='static_landmask')
C:\Users\45107\miniforge3\Lib\site-packages\numpy\lib\function_base.py:1452: RuntimeWarning: invalid value encountered in subtract
  a = op(a[slice1], a[slice2])
[DEBUG: 2023-12-14 17:50:21 : satpy.resample] Computing kd-tree parameters
[DEBUG: 2023-12-14 17:50:21 : satpy.resample] Resampling band_data
[WARNING: 2023-12-14 17:50:21 : pyresample.kd_tree] Fill value incompatible with integer data using 255 instead.
Traceback (most recent call last):
  File "C:\Users\45107\PycharmProjects\Satellite\bcdef.py", line 19, in <module>
    scn2 = scn2.resample(FLDK_area, resampler='nearest')
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 983, in resample
    new_scn.generate_possible_composites(unload)
  File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1508, in generate_possible_composites
    keepables = self._generate_composites_from_loaded_datasets()
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1527, in _generate_composites_from_loaded_datasets
    return self._generate_composites_nodes_from_loaded_datasets(needed_comp_nodes)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1533, in _generate_composites_nodes_from_loaded_datasets
    self._generate_composite(node, keepables)
  File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1591, in _generate_composite
    composite = compositor(prereq_datasets,
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\composites\__init__.py", line 1839, in __call__
    res = res.where(lsm.isin(self.values_land), res_sea)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\common.py", line 1181, in where
    return ops.where_method(self, cond, other)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\ops.py", line 178, in where_method
    return apply_ufunc(
           ^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\computation.py", line 1266, in apply_ufunc
    return apply_dataarray_vfunc(
           ^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\computation.py", line 293, in apply_dataarray_vfunc
    deep_align(
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 952, in deep_align
    aligned = align(
              ^^^^^^
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 888, in align
    aligner.align()
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 574, in align
    self.align_indexes()
  File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 421, in align_indexes
    raise ValueError(
ValueError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'y' ('y',)

Above is one of those errors. I begin to believe that @simonrp84 is right, this is some issues related to that topic. All three failures are involved with StaticImage and it happened during resampling.

I could reproduce your initial error which happened when processing several images. I have fixed that with the latest commit and I do not get the error anymore.

Wrt. to the other error it's hard to debug since you use you're own recipes. But can you double check that your land-sea mask look ok?

@yukaribbba
Copy link
Contributor

@yukaribbba Can you please also provide the code that triggers the error?

from satpy import config, Scene, find_files_and_readers, utils

from pyresample.geometry import AreaDefinition

config.set(config_path=['D:/satpy_config'])

utils.debug_on()

files = find_files_and_readers(base_dir='C:/Users/45107/Downloads/Sat/Geo/H08_20220101_1400_FLDK', reader='ahi_hsd')

FLDK_area = AreaDefinition(area_id='fldk', description='full disk', proj_id='geos',

                       projection='+proj=geos +h=35785863 +lon_0=140.7 +a=6378137.0 +rf=298.257024882273',
                       width=5500, height=5500, area_extent=(-5499999.968, -5499999.968, 5499999.968, 5499999.968))

scn1 = Scene(filenames=files)

scn1.load(['high_cloud_alpha_2000'])

scn1 = scn1.resample(FLDK_area, resampler='nearest')

scn1.show('high_cloud_alpha_2000')

scn2 = Scene(filenames=files)

scn2.load(['low_cloud_alpha_2000'])

scn2 = scn2.resample(FLDK_area, resampler='nearest')

scn2.show('low_cloud_alpha_2000')

scn3 = Scene(filenames=files)

scn3.load(['low_cloud_alpha_2000_with_background'])

scn3 = scn3.resample(FLDK_area, resampler='nearest')

scn3.show('low_cloud_alpha_2000_with_background')

scn4 = Scene(filenames=files)

scn4.load(['night_cloud_alpha_2000_with_background'])

scn4 = scn4.resample(FLDK_area, resampler='nearest')

scn4.show('night_cloud_alpha_2000_with_background')

I resample the scenes so that they can match the same area extent with landmask/blackmarble. And It's getting worse now. Only scn1 succeeded, all other three failed.

[DEBUG: 2023-12-14 17:50:19 : satpy.readers.yaml_reader] Requested orientation for Dataset B07 is 'native' (default). No flipping is applied.

[DEBUG: 2023-12-14 17:50:19 : satpy.readers.yaml_reader] Reading ('C:\Users\45107\miniforge3\Lib\site-packages\satpy\etc\readers\generic_image.yaml',)

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] PROJ data files are available at built-in paths.

[DEBUG: 2023-12-14 17:50:19 : satpy.readers.yaml_reader] Assigning to generic_image: [<FSFile "C:/Users/45107/Downloads/Sat/Geo/satpy_landmask.tif">]

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entering env context: <rasterio.env.Env object at 0x000001D553BE30D0>

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Starting outermost env

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] No GDAL environment exists

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] New GDAL environment <rasterio._env.GDALEnv object at 0x000001D5014CDF00> created

[DEBUG: 2023-12-14 17:50:19 : rasterio._filepath] Installing FilePath filesystem handler plugin...

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] GDAL_DATA found in environment.

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] PROJ data files are available at built-in paths.

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x000001D5014CDF00>.

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entered env context: <rasterio.env.Env object at 0x000001D553BE30D0>

[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Sharing flag: 0

[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Nodata success: 0, Nodata value: 0.000000

[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Dataset is started.

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting env context: <rasterio.env.Env object at 0x000001D553BE30D0>

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Cleared existing <rasterio._env.GDALEnv object at 0x000001D5014CDF00> options

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Stopped GDALEnv <rasterio._env.GDALEnv object at 0x000001D5014CDF00>.

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting outermost env

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exited env context: <rasterio.env.Env object at 0x000001D553BE30D0>

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entering env context: <rasterio.env.Env object at 0x000001D502DA3910>

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Starting outermost env

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] No GDAL environment exists

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] New GDAL environment <rasterio._env.GDALEnv object at 0x000001D502DA8520> created

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] GDAL_DATA found in environment.

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] PROJ data files are available at built-in paths.

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Started GDALEnv: self=<rasterio._env.GDALEnv object at 0x000001D502DA8520>.

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Entered env context: <rasterio.env.Env object at 0x000001D502DA3910>

[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Sharing flag: 0

[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Nodata success: 0, Nodata value: 0.000000

[DEBUG: 2023-12-14 17:50:19 : rasterio._base] Dataset is started.

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting env context: <rasterio.env.Env object at 0x000001D502DA3910>

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Cleared existing <rasterio._env.GDALEnv object at 0x000001D502DA8520> options

[DEBUG: 2023-12-14 17:50:19 : rasterio._env] Stopped GDALEnv <rasterio._env.GDALEnv object at 0x000001D502DA8520>.

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exiting outermost env

[DEBUG: 2023-12-14 17:50:19 : rasterio.env] Exited env context: <rasterio.env.Env object at 0x000001D502DA3910>

C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\dataset.py:278: UserWarning: The specified chunks separate the stored chunks along dimension "x" starting at index 4096. This could degrade performance. Instead, consider rechunking after loading.

warnings.warn(

[DEBUG: 2023-12-14 17:50:19 : satpy.composites.config_loader] Looking for composites config file images.yaml

[DEBUG: 2023-12-14 17:50:19 : satpy.composites.config_loader] No composite config found called images.yaml

[DEBUG: 2023-12-14 17:50:19 : satpy.readers.generic_image] Reading 'image.'

[WARNING: 2023-12-14 17:50:19 : satpy.readers.generic_image] cannot convert float NaN to integer

[DEBUG: 2023-12-14 17:50:19 : satpy.composites] Not all areas are the same in 'low_cloud_alpha_2000'

[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Delaying generation of DataID(name='low_cloud_alpha_2000') because of incompatible areas

[WARNING: 2023-12-14 17:50:19 : satpy.scene] The following datasets were not created and may require resampling to be generated: DataID(name='low_cloud_alpha_2000')

[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Unloading dataset: DataID(name='B07', wavelength=WavelengthRange(min=3.7, central=3.9, max=4.1, unit='µm'), resolution=2000, calibration=<2>, modifiers=())

[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Resampling DataID(name='B13', wavelength=WavelengthRange(min=10.2, central=10.4, max=10.6, unit='µm'), resolution=2000, calibration=<2>, modifiers=())

[DEBUG: 2023-12-14 17:50:19 : pyresample.geometry] Projections for data and slice areas are identical: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378137,298.257024882273,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Geostationary Satellite (Sweep Y)"],PARAMETER["Longitude of natural origin",140.7,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Satellite Height",35785863,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]]]]

[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Resampling DataID(name='_low_cloud_alpha_2000_dep_0', resolution=2000)

[DEBUG: 2023-12-14 17:50:19 : satpy.scene] Resampling DataID(name='static_landmask')

C:\Users\45107\miniforge3\Lib\site-packages\numpy\lib\function_base.py:1452: RuntimeWarning: invalid value encountered in subtract

a = op(a[slice1], a[slice2])

[DEBUG: 2023-12-14 17:50:21 : satpy.resample] Computing kd-tree parameters

[DEBUG: 2023-12-14 17:50:21 : satpy.resample] Resampling band_data

[WARNING: 2023-12-14 17:50:21 : pyresample.kd_tree] Fill value incompatible with integer data using 255 instead.

Traceback (most recent call last):

File "C:\Users\45107\PycharmProjects\Satellite\bcdef.py", line 19, in

scn2 = scn2.resample(FLDK_area, resampler='nearest')
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 983, in resample

new_scn.generate_possible_composites(unload)

File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1508, in generate_possible_composites

keepables = self._generate_composites_from_loaded_datasets()
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1527, in _generate_composites_from_loaded_datasets

return self._generate_composites_nodes_from_loaded_datasets(needed_comp_nodes)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1533, in _generate_composites_nodes_from_loaded_datasets

self._generate_composite(node, keepables)

File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\scene.py", line 1591, in _generate_composite

composite = compositor(prereq_datasets,
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\satpy\composites_init_.py", line 1839, in call

res = res.where(lsm.isin(self.values_land), res_sea)
      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\common.py", line 1181, in where

return ops.where_method(self, cond, other)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\ops.py", line 178, in where_method

return apply_ufunc(
       ^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\computation.py", line 1266, in apply_ufunc

return apply_dataarray_vfunc(
       ^^^^^^^^^^^^^^^^^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\computation.py", line 293, in apply_dataarray_vfunc

deep_align(

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 952, in deep_align

aligned = align(
          ^^^^^^

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 888, in align

aligner.align()

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 574, in align

self.align_indexes()

File "C:\Users\45107\miniforge3\Lib\site-packages\xarray\core\alignment.py", line 421, in align_indexes

raise ValueError(

ValueError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'y' ('y',)

Above is one of those errors. I begin to believe that @simonrp84 is right, this is some issues related to that topic. All three failures are involved with StaticImage and it happened during resampling.

I could reproduce your initial error which happened when processing several images. I have fixed that with the latest commit and I do not get the error anymore.

Wrt. to the other error it's hard to debug since you use you're own recipes. But can you double check that your land-sea mask look ok?

I already have produced a low-cloud before using my own landmask, which is actually from that NASA link I provided. I haven't tested your latest update and will do it tomorrow. Thanks!

@strandgren
Copy link
Collaborator Author

So I have finished the parts that were missing now and here are a few examples of the results using this PR at its current state:

ABI (GEOS-East)

ABI_GeoColor_goes_east_abi_f_2km_20231212-20231213_FR-6.mp4

AHI

AHI_GeoColor_himawari_ahi_fes_2km_20231212-20231213_FR-6.mp4

FCI

FCI_GeoColor_mtg_fci_fdss_2km_20230318-20230319_FR-6.mp4

FCI FDHSI and HRFI data resampled onto 1km grid
At reduced resolution and jpeg compressed (to fit GitHub max size limit)
FCI_GeoColor_FDHSI-HRFI_1km_reduced_and_compressed

Full resolution: https://sftp.eumetsat.int/public/file/byvy1paveeghqg20wkgvew/FCI_GeoColor_FDHSI-HRFI_1km.png

Some changes in implementation compared to imagery previously posted include:

  • Updated default configurations for night-time layers based on feedback from CIRA.
  • Update of sunz_reducion for AHI and FCI to increase true_color brightness at high sunz angels (Update Sun-zenith reducer defaults #2653).
  • Modification of sunz range for Day-Night blending in order to gain more true_color imagery.
  • Use of crude stretch for night-time layers in order to have better contrast and less saturation.

For FCI we will have to update the configurations for the low-level cloud layer when the L1C data and calibration coefficients are updated. However, I'd leave that for a future small PR.

The issue with false alarms over arid surfaces for the low-level cloud layer is due to surface emissivity differences between the two IR channels and something that we will work on to improve in the future.

@strandgren strandgren marked this pull request as ready for review December 15, 2023 07:50
Copy link
Member

@mraspaud mraspaud left a comment

Choose a reason for hiding this comment

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

LGTM!

satpy/composites/__init__.py Outdated Show resolved Hide resolved
@yukaribbba
Copy link
Contributor

@strandgren With #2690 this works pretty well. But I have a question here: I've tested the AHI data on UTC 2023-12-15 09:50, and compared my output with CIRA's. Despite the colors, there's a significant difference on LowCloud over China mainland.

My image
image

CIRA
image

Here's your land-water mask:
image

Mine:
image

@yukaribbba
Copy link
Contributor

And my composite config:

  low_cloud_alpha_2000:
    compositor: !!python/name:satpy.composites.LowCloudCompositor
    values_sea: 0
    values_land: 255
    range_land: [4.35, 6.75]
    range_sea: [1.35, 5.0]
    prerequisites:
      - compositor: !!python/name:satpy.composites.DifferenceCompositor
        prerequisites:
          - name: B13
          - name: B07
      - name: B13
      - name: static_landmask
    standard_name: low_cloud_default_spec

@strandgren
Copy link
Collaborator Author

strandgren commented Dec 15, 2023

And my composite config:

  low_cloud_alpha_2000:
    compositor: !!python/name:satpy.composites.LowCloudCompositor
    values_sea: 0
    values_land: 255
    range_land: [4.35, 6.75]
    range_sea: [1.35, 5.0]
    prerequisites:
      - compositor: !!python/name:satpy.composites.DifferenceCompositor
        prerequisites:
          - name: B13
          - name: B07
      - name: B13
      - name: static_landmask
    standard_name: low_cloud_default_spec

remove range_land and range_sea here. Then you will use the default ones which are the ones provided by CIRA. That should hopefully yield more similar results. The values above we only use for FCI at the moment during commissioning where the channels are not yet properly calibrated.

Please note that in one of the recent commits I changes the term sea to water so if you want to define your own range for water surface types, you'd have to use range_water instead.

@yukaribbba
Copy link
Contributor

And my composite config:

  low_cloud_alpha_2000:
    compositor: !!python/name:satpy.composites.LowCloudCompositor
    values_sea: 0
    values_land: 255
    range_land: [4.35, 6.75]
    range_sea: [1.35, 5.0]
    prerequisites:
      - compositor: !!python/name:satpy.composites.DifferenceCompositor
        prerequisites:
          - name: B13
          - name: B07
      - name: B13
      - name: static_landmask
    standard_name: low_cloud_default_spec

remove range_land and range_sea here. Then you will use the default ones which are the ones provided by CIRA. That should hopefully yield more similar results. The values above we only use for FCI at the moment during commissioning where the channels are not yet properly calibrated.

Please note that in one of the recent commits I changes the term sea to water so if you want to define your own range for water surface types, you'd have to use range_water instead.

Great!
image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
component:compositors enhancement code enhancements, features, improvements
Projects
Development

Successfully merging this pull request may close these issues.

6 participants