methane_super_emitters.utils

View Source

  1import glob
  2import os
  3import random
  4import shutil
  5import numpy as np
  6import netCDF4
  7import datetime
  8import geedim
  9import logging
 10import pathlib
 11import rasterio
 12from rasterio.plot import reshape_as_image
 13from rasterio.transform import from_origin
 14from PIL import Image
 15from dateutil.relativedelta import relativedelta
 16
 17
 18def s2_download(npz_file, output_dir):
 19    """Download a Sentinel-2 image corresponding to a TROPOMI patch.
 20
 21    Parameters
 22    ----------
 23    npz_file: str
 24        A patch from a TROPOMI file used for training or inference in our model.
 25    output_dir: str
 26        Output directory to store the Sentinel-2 data.
 27    """
 28    geedim.Initialize()
 29    data = np.load(npz_file, allow_pickle=True)
 30    lat = data["lat_bounds"]
 31    lon = data["lon_bounds"]
 32    lat_min = float(lat.min())
 33    lat_max = float(lat.max())
 34    lon_min = float(lon.min())
 35    lon_max = float(lon.max())
 36    polygon = {
 37        "type": "Polygon",
 38        "coordinates": [
 39            [
 40                [lat_min, lon_min],
 41                [lat_min, lon_max],
 42                [lat_max, lon_max],
 43                [lat_max, lon_min],
 44                [lat_min, lon_min],
 45            ]
 46        ],
 47    }
 48    coll = geedim.MaskedCollection.from_name("COPERNICUS/S2_HARMONIZED")
 49    date = min(data["time"])
 50    start_date = date - relativedelta(weeks=2)
 51    end_date = date + relativedelta(weeks=2)
 52    try:
 53        coll = coll.search(
 54            start_date=start_date,
 55            end_date=end_date,
 56            region=polygon,
 57            cloudless_portion=0.5,
 58        )
 59        comp_im = coll.composite(method="mosaic", region=polygon)
 60        stem = pathlib.Path(npz_file).stem
 61        output_file = os.path.join(output_dir, f"{stem}.tif")
 62        comp_im.download(
 63            output_file,
 64            region=polygon,
 65            crs="EPSG:4326",
 66            scale=20,
 67            overwrite=True,
 68            bands=["B4", "B3", "B2"],
 69        )
 70    except ValueError:
 71        logging.warning(
 72            f"No S2 data find matching the query constructed from {npz_file}"
 73        )
 74
 75
 76def s2_preview_png(tiff_file, output_dir):
 77    """Create a PNG preview of a S2 GeoTIFF file.
 78
 79    Parameters
 80    ----------
 81    tiff_file: str
 82        A GeoTIFF with the Sentinel-2 scene.
 83    output_dir: str
 84        Output directory for the result.
 85    """
 86    stem = pathlib.Path(tiff_file).stem
 87    with rasterio.open(tiff_file, "r") as src:
 88        rgb = src.read([1, 2, 3])
 89        rgb = (255 * (rgb / np.max(rgb))).astype(np.uint8)
 90        rgb_image = reshape_as_image(rgb)
 91        Image.fromarray(rgb_image).save(os.path.join(output_dir, f"{stem}.png"))
 92
 93
 94def npz_to_geotiff(npz_file, output_dir):
 95    """Store methane concentrations in a GeoTIFF for preview.
 96
 97    Parameters
 98    ----------
 99    npz_file: str
100        An NPZ file as used for training and inference.
101    output_dir: str
102        Output directory for the result.
103    """
104    stem = pathlib.Path(npz_file).stem
105    data = np.load(npz_file)
106    ch4 = data["methane"]
107    mask = data["mask"]
108    lat = data["lat_bounds"]
109    lon = data["lon_bounds"]
110    lat_res = (lat.max() - lat.min()) / lat.shape[0]
111    lon_res = (lon.max() - lon.min()) / lat.shape[1]
112    transform = from_origin(lon.min(), lat.max(), lon_res, lat_res)
113    ch4[mask] = -9999
114    with rasterio.open(
115        os.path.join(output_dir, f"{stem}_ch4.tif"),
116        "w",
117        driver="GTiff",
118        height=ch4.shape[0],
119        width=ch4.shape[1],
120        count=1,
121        dtype=np.float32,
122        crs="EPSG:4326",
123        transform=transform,
124        nodata=-9999,
125    ) as dst:
126        dst.write(ch4, 1)
127
128
129def sample_files(glob_pattern, output_dir, n):
130    """A small utility function to sample files from a directory at random.
131
132    Parameters
133    ----------
134    glob_pattern: str
135        A string of the form /path/to/*.npz or so. To select files from.
136    output_dir: str
137        A directory name into which to copy files.
138    n: int
139        How many files to take at random.
140    """
141    for path in random.sample(list(glob.glob(glob_pattern)), k=n):
142        shutil.copy(path, output_dir)
143
144
145def destripe(fd):
146    """A destriping algorithm for TROPOMI data.
147
148    Parameters
149    ----------
150    fd: a NetCDF4 file handler
151
152    Returns
153    -------
154    A numpy array with destriped methane ppm amounts.
155    """
156    ch4 = fd["/PRODUCT/methane_mixing_ratio"][:]
157    ch4corr = fd["/PRODUCT/methane_mixing_ratio_bias_corrected"][:]
158    ch4corr[ch4corr > 1e20] = np.nan
159    ch4corrdestrip = ch4corr.copy() * np.nan
160    # get the number of rows
161    n = ch4corr.shape[1]
162    # get the number of columns
163    m = ch4corr.shape[2]
164    back = np.zeros((n, m)) * np.nan
165    for i in range(m):
166        # define half window size
167        ws = 7
168        if i < ws:
169            st = 0
170            sp = i + ws
171        elif m - i < ws:
172            st = i - ws
173            sp = m - 1
174        else:
175            st = i - ws
176            sp = i + ws
177        back[:, i] = np.nanmedian(ch4corr[0, :, st:sp], axis=1)
178    this = ch4corr[0, :, :] - back
179    stripes = np.zeros((n, m)) * np.nan
180    for j in range(n):
181        ws = 60
182        if j < ws:
183            st = 0
184            sp = j + ws
185        elif n - j < ws:
186            st = j - ws
187            sp = n - 1
188        else:
189            st = j - ws
190            sp = j + ws
191        stripes[j, :] = np.nanmedian(this[st:sp, :], axis=0)
192    ch4corrdestrip[0, :, :] = ch4corr[0, :, :] - stripes
193    return ch4corrdestrip
194
195
196def parse_date(date_str, time_str):
197    """Parse date from a simple string format.
198
199    Parameters
200    ----------
201    date_str: str
202        Date of the form YYYYMMDD
203    time_str: str
204        Time of the form HH:MM:SS
205
206    Returns
207    -------
208    datetime.datetime
209    """
210    return datetime.datetime.strptime(date_str + time_str, "%Y%m%d%H:%M:%S")
211
212
213def patch_generator(file_path):
214    """A generator for 32x32 pixel patches with 50% overlap.
215
216    Parameters
217    ----------
218    file_path: str
219        A path to a TROPOMI netCDF4 file.
220
221    Yields
222    ------
223    dict
224        A dictionary with patch arrays.
225    """
226    with netCDF4.Dataset(file_path, "r") as fd:
227        destriped = destripe(fd)
228        rows = destriped.shape[1]
229        cols = destriped.shape[2]
230        for row in range(0, rows, 16):
231            for col in range(0, cols, 16):
232                if row + 32 < rows and col + 32 < cols:
233                    methane_window = destriped[0][row : row + 32][:, col : col + 32]
234                    original = fd["PRODUCT/methane_mixing_ratio_bias_corrected"][:][0][
235                        row : row + 32
236                    ][:, col : col + 32]
237                    lat_window = fd["PRODUCT/latitude"][:][0][row : row + 32][
238                        :, col : col + 32
239                    ]
240                    lon_window = fd["PRODUCT/longitude"][:][0][row : row + 32][
241                        :, col : col + 32
242                    ]
243                    lat_bounds_window = fd[
244                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds"
245                    ][:][0][row : row + 32][:, col : col + 32]
246                    lon_bounds_window = fd[
247                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds"
248                    ][:][0][row : row + 32][:, col : col + 32]
249                    qa_window = fd["PRODUCT/qa_value"][:][0][row : row + 32][
250                        :, col : col + 32
251                    ]
252                    u10_window = fd["PRODUCT/SUPPORT_DATA/INPUT_DATA/eastward_wind"][:][
253                        0
254                    ][row : row + 32][:, col : col + 32]
255                    v10_window = fd["PRODUCT/SUPPORT_DATA/INPUT_DATA/northward_wind"][
256                        :
257                    ][0][row : row + 32][:, col : col + 32]
258                    sza_window = fd[
259                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/solar_zenith_angle"
260                    ][:][0][row : row + 32][:, col : col + 32]
261                    vza_window = fd[
262                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/viewing_zenith_angle"
263                    ][:][0][row : row + 32][:, col : col + 32]
264                    saa_window = fd[
265                        "/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/solar_azimuth_angle"
266                    ][:][0][row : row + 32][:, col : col + 32]
267                    vaa_window = fd[
268                        "/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/viewing_azimuth_angle"
269                    ][:][0][row : row + 32][:, col : col + 32]
270                    f_delta = saa_window - vaa_window
271                    raa_window = np.where(
272                        f_delta < 0, np.abs(f_delta + 180.0), np.abs(f_delta - 180.0)
273                    )
274                    raa_window = abs(180.0 - raa_window)
275                    scattering_angle = -np.cos(np.radians(sza_window)) * np.cos(
276                        np.radians(vza_window)
277                    ) + np.sin(np.radians(sza_window)) * np.sin(
278                        np.radians(vza_window)
279                    ) * np.cos(
280                        np.radians(raa_window)
281                    )
282                    sa_std_window = fd[
283                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_altitude_precision"
284                    ][:][0][row : row + 32][:, col : col + 32]
285                    cloud_fraction_window = fd[
286                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/cloud_fraction_VIIRS_SWIR_IFOV"
287                    ][:][0][row : row + 32][:, col : col + 32]
288                    cirrus_reflectance_window = fd[
289                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/reflectance_cirrus_VIIRS_SWIR"
290                    ][:][0][row : row + 32][:, col : col + 32]
291                    methane_ratio_std_window = fd[
292                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/methane_ratio_weak_strong_standard_deviation"
293                    ][:][0][row : row + 32][:, col : col + 32]
294                    methane_precision_window = fd[
295                        "PRODUCT/methane_mixing_ratio_precision"
296                    ][:][0][row : row + 32][:, col : col + 32]
297                    surface_albedo_window = fd[
298                        "PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/surface_albedo_SWIR"
299                    ][:][0][row : row + 32][:, col : col + 32]
300                    surface_albedo_precision_window = fd[
301                        "PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/surface_albedo_SWIR_precision"
302                    ][:][0][row : row + 32][:, col : col + 32]
303                    aerosol_optical_thickness_window = fd[
304                        "PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/aerosol_optical_thickness_SWIR"
305                    ][:][0][row : row + 32][:, col : col + 32]
306                    times = fd["PRODUCT/time_utc"][:][0][row : row + 32]
307                    try:
308                        parsed_time = [
309                            datetime.datetime.fromisoformat(time_str[:19])
310                            for time_str in times
311                        ]
312                    except ValueError:
313                        continue
314                    if original.mask.sum() < 0.2 * 32 * 32:
315                        yield {
316                            "methane": methane_window,
317                            "lat": lat_window,
318                            "lon": lon_window,
319                            "lat_bounds": lat_bounds_window,
320                            "lon_bounds": lon_bounds_window,
321                            "qa": qa_window,
322                            "time": parsed_time,
323                            "mask": original.mask,
324                            "non_destriped": original,
325                            "u10": u10_window,
326                            "v10": v10_window,
327                            "sza": np.cos(np.radians(sza_window)),
328                            "vza": np.cos(np.radians(vza_window)),
329                            "scattering_angle": scattering_angle,
330                            "sa_std": sa_std_window,
331                            "cloud_fraction": cloud_fraction_window,
332                            "cirrus_reflectance": cirrus_reflectance_window,
333                            "methane_ratio_std": methane_ratio_std_window,
334                            "methane_precision": methane_precision_window,
335                            "surface_albedo": surface_albedo_window,
336                            "surface_albedo_precision": surface_albedo_precision_window,
337                            "aerosol_optical_thickness": aerosol_optical_thickness_window,
338                        }

def s2_download(npz_file, output_dir): View Source

19def s2_download(npz_file, output_dir):
20    """Download a Sentinel-2 image corresponding to a TROPOMI patch.
21
22    Parameters
23    ----------
24    npz_file: str
25        A patch from a TROPOMI file used for training or inference in our model.
26    output_dir: str
27        Output directory to store the Sentinel-2 data.
28    """
29    geedim.Initialize()
30    data = np.load(npz_file, allow_pickle=True)
31    lat = data["lat_bounds"]
32    lon = data["lon_bounds"]
33    lat_min = float(lat.min())
34    lat_max = float(lat.max())
35    lon_min = float(lon.min())
36    lon_max = float(lon.max())
37    polygon = {
38        "type": "Polygon",
39        "coordinates": [
40            [
41                [lat_min, lon_min],
42                [lat_min, lon_max],
43                [lat_max, lon_max],
44                [lat_max, lon_min],
45                [lat_min, lon_min],
46            ]
47        ],
48    }
49    coll = geedim.MaskedCollection.from_name("COPERNICUS/S2_HARMONIZED")
50    date = min(data["time"])
51    start_date = date - relativedelta(weeks=2)
52    end_date = date + relativedelta(weeks=2)
53    try:
54        coll = coll.search(
55            start_date=start_date,
56            end_date=end_date,
57            region=polygon,
58            cloudless_portion=0.5,
59        )
60        comp_im = coll.composite(method="mosaic", region=polygon)
61        stem = pathlib.Path(npz_file).stem
62        output_file = os.path.join(output_dir, f"{stem}.tif")
63        comp_im.download(
64            output_file,
65            region=polygon,
66            crs="EPSG:4326",
67            scale=20,
68            overwrite=True,
69            bands=["B4", "B3", "B2"],
70        )
71    except ValueError:
72        logging.warning(
73            f"No S2 data find matching the query constructed from {npz_file}"
74        )

Download a Sentinel-2 image corresponding to a TROPOMI patch.

Parameters

npz_file (str): A patch from a TROPOMI file used for training or inference in our model.
output_dir (str): Output directory to store the Sentinel-2 data.

def s2_preview_png(tiff_file, output_dir): View Source

77def s2_preview_png(tiff_file, output_dir):
78    """Create a PNG preview of a S2 GeoTIFF file.
79
80    Parameters
81    ----------
82    tiff_file: str
83        A GeoTIFF with the Sentinel-2 scene.
84    output_dir: str
85        Output directory for the result.
86    """
87    stem = pathlib.Path(tiff_file).stem
88    with rasterio.open(tiff_file, "r") as src:
89        rgb = src.read([1, 2, 3])
90        rgb = (255 * (rgb / np.max(rgb))).astype(np.uint8)
91        rgb_image = reshape_as_image(rgb)
92        Image.fromarray(rgb_image).save(os.path.join(output_dir, f"{stem}.png"))

Create a PNG preview of a S2 GeoTIFF file.

Parameters

tiff_file (str): A GeoTIFF with the Sentinel-2 scene.
output_dir (str): Output directory for the result.

def npz_to_geotiff(npz_file, output_dir): View Source

 95def npz_to_geotiff(npz_file, output_dir):
 96    """Store methane concentrations in a GeoTIFF for preview.
 97
 98    Parameters
 99    ----------
100    npz_file: str
101        An NPZ file as used for training and inference.
102    output_dir: str
103        Output directory for the result.
104    """
105    stem = pathlib.Path(npz_file).stem
106    data = np.load(npz_file)
107    ch4 = data["methane"]
108    mask = data["mask"]
109    lat = data["lat_bounds"]
110    lon = data["lon_bounds"]
111    lat_res = (lat.max() - lat.min()) / lat.shape[0]
112    lon_res = (lon.max() - lon.min()) / lat.shape[1]
113    transform = from_origin(lon.min(), lat.max(), lon_res, lat_res)
114    ch4[mask] = -9999
115    with rasterio.open(
116        os.path.join(output_dir, f"{stem}_ch4.tif"),
117        "w",
118        driver="GTiff",
119        height=ch4.shape[0],
120        width=ch4.shape[1],
121        count=1,
122        dtype=np.float32,
123        crs="EPSG:4326",
124        transform=transform,
125        nodata=-9999,
126    ) as dst:
127        dst.write(ch4, 1)

Store methane concentrations in a GeoTIFF for preview.

Parameters

npz_file (str): An NPZ file as used for training and inference.
output_dir (str): Output directory for the result.

def sample_files(glob_pattern, output_dir, n): View Source

130def sample_files(glob_pattern, output_dir, n):
131    """A small utility function to sample files from a directory at random.
132
133    Parameters
134    ----------
135    glob_pattern: str
136        A string of the form /path/to/*.npz or so. To select files from.
137    output_dir: str
138        A directory name into which to copy files.
139    n: int
140        How many files to take at random.
141    """
142    for path in random.sample(list(glob.glob(glob_pattern)), k=n):
143        shutil.copy(path, output_dir)

A small utility function to sample files from a directory at random.

Parameters

glob_pattern (str): A string of the form /path/to/*.npz or so. To select files from.
output_dir (str): A directory name into which to copy files.
n (int): How many files to take at random.

def destripe(fd): View Source

146def destripe(fd):
147    """A destriping algorithm for TROPOMI data.
148
149    Parameters
150    ----------
151    fd: a NetCDF4 file handler
152
153    Returns
154    -------
155    A numpy array with destriped methane ppm amounts.
156    """
157    ch4 = fd["/PRODUCT/methane_mixing_ratio"][:]
158    ch4corr = fd["/PRODUCT/methane_mixing_ratio_bias_corrected"][:]
159    ch4corr[ch4corr > 1e20] = np.nan
160    ch4corrdestrip = ch4corr.copy() * np.nan
161    # get the number of rows
162    n = ch4corr.shape[1]
163    # get the number of columns
164    m = ch4corr.shape[2]
165    back = np.zeros((n, m)) * np.nan
166    for i in range(m):
167        # define half window size
168        ws = 7
169        if i < ws:
170            st = 0
171            sp = i + ws
172        elif m - i < ws:
173            st = i - ws
174            sp = m - 1
175        else:
176            st = i - ws
177            sp = i + ws
178        back[:, i] = np.nanmedian(ch4corr[0, :, st:sp], axis=1)
179    this = ch4corr[0, :, :] - back
180    stripes = np.zeros((n, m)) * np.nan
181    for j in range(n):
182        ws = 60
183        if j < ws:
184            st = 0
185            sp = j + ws
186        elif n - j < ws:
187            st = j - ws
188            sp = n - 1
189        else:
190            st = j - ws
191            sp = j + ws
192        stripes[j, :] = np.nanmedian(this[st:sp, :], axis=0)
193    ch4corrdestrip[0, :, :] = ch4corr[0, :, :] - stripes
194    return ch4corrdestrip

A destriping algorithm for TROPOMI data.

Parameters

fd (a NetCDF4 file handler):

Returns

A numpy array with destriped methane ppm amounts.

def parse_date(date_str, time_str): View Source

197def parse_date(date_str, time_str):
198    """Parse date from a simple string format.
199
200    Parameters
201    ----------
202    date_str: str
203        Date of the form YYYYMMDD
204    time_str: str
205        Time of the form HH:MM:SS
206
207    Returns
208    -------
209    datetime.datetime
210    """
211    return datetime.datetime.strptime(date_str + time_str, "%Y%m%d%H:%M:%S")

Parse date from a simple string format.

Parameters

date_str (str): Date of the form YYYYMMDD
time_str (str): Time of the form HH:MM:SS

Returns

datetime.datetime

def patch_generator(file_path): View Source

214def patch_generator(file_path):
215    """A generator for 32x32 pixel patches with 50% overlap.
216
217    Parameters
218    ----------
219    file_path: str
220        A path to a TROPOMI netCDF4 file.
221
222    Yields
223    ------
224    dict
225        A dictionary with patch arrays.
226    """
227    with netCDF4.Dataset(file_path, "r") as fd:
228        destriped = destripe(fd)
229        rows = destriped.shape[1]
230        cols = destriped.shape[2]
231        for row in range(0, rows, 16):
232            for col in range(0, cols, 16):
233                if row + 32 < rows and col + 32 < cols:
234                    methane_window = destriped[0][row : row + 32][:, col : col + 32]
235                    original = fd["PRODUCT/methane_mixing_ratio_bias_corrected"][:][0][
236                        row : row + 32
237                    ][:, col : col + 32]
238                    lat_window = fd["PRODUCT/latitude"][:][0][row : row + 32][
239                        :, col : col + 32
240                    ]
241                    lon_window = fd["PRODUCT/longitude"][:][0][row : row + 32][
242                        :, col : col + 32
243                    ]
244                    lat_bounds_window = fd[
245                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds"
246                    ][:][0][row : row + 32][:, col : col + 32]
247                    lon_bounds_window = fd[
248                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds"
249                    ][:][0][row : row + 32][:, col : col + 32]
250                    qa_window = fd["PRODUCT/qa_value"][:][0][row : row + 32][
251                        :, col : col + 32
252                    ]
253                    u10_window = fd["PRODUCT/SUPPORT_DATA/INPUT_DATA/eastward_wind"][:][
254                        0
255                    ][row : row + 32][:, col : col + 32]
256                    v10_window = fd["PRODUCT/SUPPORT_DATA/INPUT_DATA/northward_wind"][
257                        :
258                    ][0][row : row + 32][:, col : col + 32]
259                    sza_window = fd[
260                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/solar_zenith_angle"
261                    ][:][0][row : row + 32][:, col : col + 32]
262                    vza_window = fd[
263                        "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/viewing_zenith_angle"
264                    ][:][0][row : row + 32][:, col : col + 32]
265                    saa_window = fd[
266                        "/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/solar_azimuth_angle"
267                    ][:][0][row : row + 32][:, col : col + 32]
268                    vaa_window = fd[
269                        "/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/viewing_azimuth_angle"
270                    ][:][0][row : row + 32][:, col : col + 32]
271                    f_delta = saa_window - vaa_window
272                    raa_window = np.where(
273                        f_delta < 0, np.abs(f_delta + 180.0), np.abs(f_delta - 180.0)
274                    )
275                    raa_window = abs(180.0 - raa_window)
276                    scattering_angle = -np.cos(np.radians(sza_window)) * np.cos(
277                        np.radians(vza_window)
278                    ) + np.sin(np.radians(sza_window)) * np.sin(
279                        np.radians(vza_window)
280                    ) * np.cos(
281                        np.radians(raa_window)
282                    )
283                    sa_std_window = fd[
284                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_altitude_precision"
285                    ][:][0][row : row + 32][:, col : col + 32]
286                    cloud_fraction_window = fd[
287                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/cloud_fraction_VIIRS_SWIR_IFOV"
288                    ][:][0][row : row + 32][:, col : col + 32]
289                    cirrus_reflectance_window = fd[
290                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/reflectance_cirrus_VIIRS_SWIR"
291                    ][:][0][row : row + 32][:, col : col + 32]
292                    methane_ratio_std_window = fd[
293                        "PRODUCT/SUPPORT_DATA/INPUT_DATA/methane_ratio_weak_strong_standard_deviation"
294                    ][:][0][row : row + 32][:, col : col + 32]
295                    methane_precision_window = fd[
296                        "PRODUCT/methane_mixing_ratio_precision"
297                    ][:][0][row : row + 32][:, col : col + 32]
298                    surface_albedo_window = fd[
299                        "PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/surface_albedo_SWIR"
300                    ][:][0][row : row + 32][:, col : col + 32]
301                    surface_albedo_precision_window = fd[
302                        "PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/surface_albedo_SWIR_precision"
303                    ][:][0][row : row + 32][:, col : col + 32]
304                    aerosol_optical_thickness_window = fd[
305                        "PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/aerosol_optical_thickness_SWIR"
306                    ][:][0][row : row + 32][:, col : col + 32]
307                    times = fd["PRODUCT/time_utc"][:][0][row : row + 32]
308                    try:
309                        parsed_time = [
310                            datetime.datetime.fromisoformat(time_str[:19])
311                            for time_str in times
312                        ]
313                    except ValueError:
314                        continue
315                    if original.mask.sum() < 0.2 * 32 * 32:
316                        yield {
317                            "methane": methane_window,
318                            "lat": lat_window,
319                            "lon": lon_window,
320                            "lat_bounds": lat_bounds_window,
321                            "lon_bounds": lon_bounds_window,
322                            "qa": qa_window,
323                            "time": parsed_time,
324                            "mask": original.mask,
325                            "non_destriped": original,
326                            "u10": u10_window,
327                            "v10": v10_window,
328                            "sza": np.cos(np.radians(sza_window)),
329                            "vza": np.cos(np.radians(vza_window)),
330                            "scattering_angle": scattering_angle,
331                            "sa_std": sa_std_window,
332                            "cloud_fraction": cloud_fraction_window,
333                            "cirrus_reflectance": cirrus_reflectance_window,
334                            "methane_ratio_std": methane_ratio_std_window,
335                            "methane_precision": methane_precision_window,
336                            "surface_albedo": surface_albedo_window,
337                            "surface_albedo_precision": surface_albedo_precision_window,
338                            "aerosol_optical_thickness": aerosol_optical_thickness_window,
339                        }

A generator for 32x32 pixel patches with 50% overlap.

Parameters

file_path (str): A path to a TROPOMI netCDF4 file.

Yields

dict: A dictionary with patch arrays.