# Plot selected MODIS cloud-property overpasses over SGP # Panels: LWP, CDNC, COT, CER, CTH # For each MODIS file, this writes two figures: # 1) no pixel-quality filters applied, except fill-value removal # 2) updated SGP-style local/SZA/VZA/phase/SPI/CER/COT filters # + final sampling constraints applied import os from glob import glob import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from matplotlib.colors import Normalize from matplotlib.axes import Axes import cartopy import cartopy.crs as ccrs from cartopy.mpl.geoaxes import GeoAxes from netCDF4 import Dataset from pyhdf import SD from scipy.ndimage import zoom from astropy.convolution import convolve, Ring2DKernel GeoAxes._pcolormesh_patched = Axes.pcolormesh # ----------------------------- # User settings # ----------------------------- lat_site=36.63 ## SGP lat lon ## lon_site=-97.47 half_width_deg = 0.50 # 0.50 means a 1° x 1° box around SGP local_filter_width = 2 # updated code uses width=2 # Updated base filter thresholds from the latest uploaded code SZA_MAX = 65 VZA_MAX = 55 SPI_MAX = 30 CER_MIN = 4 COT_MIN = 4 SINGLE_LAYER_FLAG = 1 LIQUID_PHASE_IR = 1 # Additional final sampling constraints used in the combined dataframe analysis COT_MAX = 70 CF_MIN = 0.9 CER_MAX = 30 CTT_MIN = 268 LWP_MIN = 50 LWP_MAX = 300 alpha = 1.37e-6 # ----------------------------- # Helper functions # ----------------------------- def LWP_given_COT_CER(COT, CER): return 5.0 / 9.0 * 0.8 * COT * CER def CDNC_given_COT_CER(COT, CER): return alpha * COT**0.5 / (CER / 10000.0)**2.5 def loadVarFromHDF(filename, varname): fileID = SD.SD(filename) sds = fileID.select(varname) variable = np.squeeze(sds.get()) variable = np.flipud(variable) fileID.end() return variable def computeLocalFilterRemoveMask(MOD06filename, width=2, return_components=False): """ Updated local filter from the latest uploaded code. Return True for pixels to remove because a neighbouring pixel is: - an unsuccessful liquid retrieval, excluding pixels where PCL succeeded, or - ice cloud. The latest code does not include clear-sky-restoral pixels in the local-filter source mask, and it uses width=2. """ ds = Dataset(MOD06filename, "r") QA1km = (ds.variables["Quality_Assurance_1km"][:]).data[:, :-4, :].astype(np.uint8) phase1km = (QA1km[:, :, 2] & 0b00000111) >> 0 liquidPhaseCloud1km = phase1km == 2 successfulRetrieval1km = ((QA1km[:, :, 2] & 0b00001000) >> 3) == 1 # PCL success flag: do not count a liquid retrieval as unsuccessful if PCL succeeded. successfulLiquidRetrievalPCL = ((QA1km[:, :, 8] & 0b11110000) >> 4) == 10 unsuccessfulLiquidRetrieval1km = np.logical_and( ~(successfulRetrieval1km | successfulLiquidRetrievalPCL), liquidPhaseCloud1km, ) iceCloud1km = phase1km == 3 bad_neighbour_source = np.logical_or.reduce((unsuccessfulLiquidRetrieval1km, iceCloud1km)) kernel = Ring2DKernel(radius_in=0, width=width) localFilterRemove = ( convolve(bad_neighbour_source.astype(float), kernel, boundary="fill", fill_value=0.0) > 0.0 ) ds.close() if return_components: return localFilterRemove.astype(bool), unsuccessfulLiquidRetrieval1km, iceCloud1km return localFilterRemove.astype(bool) def trim_last4_columns(arr): """Trim final 4 columns so 1 km science arrays match the QA/local-filter mask.""" return arr[:, :-4] def crop_to_shape(arr, shape): """Crop an array to a target (ny, nx) shape.""" return arr[:shape[0], :shape[1]] def common_2d_shape(*arrays): """Return the smallest common 2D shape across all given arrays.""" ny = min(arr.shape[0] for arr in arrays) nx = min(arr.shape[1] for arr in arrays) return ny, nx def clean_fill(arr, fill_limit=-99.0): arr = arr.astype(float) arr[arr < fill_limit] = np.nan return arr def add_hyytiala_box(ax): ax.plot([lon_site - half_width_deg, lon_site - half_width_deg], [lat_site - half_width_deg, lat_site + half_width_deg], linestyle="--", c="#8b3300", alpha=0.925, zorder=800, linewidth=1.2, transform=ccrs.PlateCarree()) ax.plot([lon_site + half_width_deg, lon_site + half_width_deg], [lat_site - half_width_deg, lat_site + half_width_deg], linestyle="--", c="#8b3300", alpha=0.925, zorder=800, linewidth=1.2, transform=ccrs.PlateCarree()) ax.plot([lon_site - half_width_deg, lon_site + half_width_deg], [lat_site - half_width_deg, lat_site - half_width_deg], linestyle="--", c="#8b3300", alpha=0.925, zorder=800, linewidth=1.2, transform=ccrs.PlateCarree()) ax.plot([lon_site - half_width_deg, lon_site + half_width_deg], [lat_site + half_width_deg, lat_site + half_width_deg], linestyle="--", c="#8b3300", alpha=0.925, zorder=800, linewidth=1.2, transform=ccrs.PlateCarree()) ax.plot(lon_site, lat_site, marker="^", color="red", markersize=5, transform=ccrs.PlateCarree(), zorder=900) def add_base_features(ax): """Add simple land/ocean/coastline/border features to a map axis.""" ocean = cartopy.feature.NaturalEarthFeature(category="physical", name="ocean", scale="50m") coastline = cartopy.feature.NaturalEarthFeature(category="physical", name="coastline", scale="50m") land = cartopy.feature.NaturalEarthFeature(category="physical", name="land", scale="50m") borders = cartopy.feature.NaturalEarthFeature(category="cultural", name="admin_0_boundary_lines_land", scale="50m") ax.add_feature(ocean, facecolor="#dedede", edgecolor="#303030", linewidth=0.5, alpha=0.5, zorder=1) ax.add_feature(land, facecolor="#eeeeee", edgecolor="none", linewidth=0, zorder=2) ax.add_feature(coastline, linestyle="-", facecolor="none", edgecolor="#303030", linewidth=0.5, alpha=0.5, zorder=12) ax.add_feature(borders, linestyle="-", facecolor="none", edgecolor="#303030", linewidth=0.25, alpha=0.5, zorder=13) def setup_large_swath_axis(ax): """Large map for the first LWP panel, like the original full-overpass figure.""" ax.set_xlim([-1500000, 1500000]) ax.set_ylim([-1500000, 1500000]) add_base_features(ax) for ii in range(-20, 20, 2): ax.plot([int(lon_site) + ii, int(lon_site) + ii], [-89, 89], linestyle="--", c="#202020", alpha=0.25, zorder=30, linewidth=0.5, transform=ccrs.PlateCarree()) ax.plot([-170, 170], [int(lat_site) + ii, int(lat_site) + ii], linestyle="--", c="#202020", alpha=0.25, zorder=30, linewidth=0.5, transform=ccrs.PlateCarree()) add_hyytiala_box(ax) def setup_box_axis(ax): """ Rectangular local map for CDNC/COT/CER/CTH, similar to the earlier figure. The dashed brown rectangle marks the true 1° x 1° latitude-longitude sampling domain. The plotted data are not masked to that domain; the map limits only crop the display to a slightly larger local rectangular area around SGP. """ # Rectangular projected-metre view around SGP. # East-west extent is smaller than north-south extent because the 1° x 1° # lat-lon sampling box is physically tall at 61.83°N. local_half_x_m = 65000 local_half_y_m = 85000 ax.set_xlim([-local_half_x_m, local_half_x_m]) ax.set_ylim([-local_half_y_m, local_half_y_m]) add_base_features(ax) add_hyytiala_box(ax) def mask_outside_hyytiala_box(values, lat_1km, lon_1km): """Keep values only inside the 1° x 1° SGP lat-lon box.""" domain_mask = ( (lat_1km >= lat_site - half_width_deg) & (lat_1km <= lat_site + half_width_deg) & (lon_1km >= lon_site - half_width_deg) & (lon_1km <= lon_site + half_width_deg) ) out = values.copy() out[~domain_mask] = np.nan return out def plot_5panel_figure(datafile, lon_1km, lat_1km, fields, output_file, title_suffix): """ Five-panel figure: - LWP is shown on the large full-overpass/swath map. - CDNC, COT, CER and CTH are shown on rectangular local map frames around SGP. The dashed 1° x 1° sampling box is overlaid for reference, but the data are NOT masked to that box; data are shown across the whole local rectangular map panel wherever available. """ panels = [ ("LWP", fields["LWP"], "Liquid Water Path (g m$^{-2}$)", plt.get_cmap("viridis"), Normalize(0, 200), "swath"), ("CDNC", fields["CDNC"], "CDNC (cm$^{-3}$)", plt.get_cmap("PiYG"), Normalize(0, 500), "box"), ("COT", fields["COT"], "Cloud Optical Thickness", plt.get_cmap("YlGnBu"), Normalize(0, 70), "box"), ("CER", fields["CER"], "Cloud Effective Radius (µm)", plt.get_cmap("cividis"), Normalize(0, 30), "box"), ("CTH", fields["CTH"], "Cloud Top Height (km)", plt.get_cmap("magma"), Normalize(0, 10), "box"), ] fig = plt.figure(figsize=(10.5, 12.0)) # Manual positions: # - LWP keeps the large swath/overpass view. # - Other panels use rectangular local axes, closer to your earlier layout. axis_specs = [ [0.055, 0.675, 0.39, 0.255], # LWP, large swath [0.545, 0.675, 0.34, 0.255], # CDNC, rectangular local panel [0.105, 0.385, 0.34, 0.255], # COT, rectangular local panel [0.545, 0.385, 0.34, 0.255], # CER, rectangular local panel [0.105, 0.095, 0.34, 0.255], # CTH, rectangular local panel ] for idx, (panel_title, values, cbar_label, cmap, norm, view_type) in enumerate(panels): if view_type == "swath": proj = ccrs.TransverseMercator( central_longitude=float(np.nanmean(lon_1km)), central_latitude=float(np.nanmean(lat_1km)), ) else: proj = ccrs.TransverseMercator(central_longitude=lon_site, central_latitude=lat_site) ax = fig.add_axes(axis_specs[idx], projection=proj) if view_type == "swath": setup_large_swath_axis(ax) plot_values = values else: setup_box_axis(ax) # Do not mask to the 1° x 1° sampling box here. # The axis limits clip the plot to the rectangular local frame, while # the dashed rectangle only marks the sampling domain. plot_values = values im = ax.pcolormesh( lon_1km, lat_1km, plot_values, cmap=cmap, norm=norm, shading="nearest", zorder=50, transform=ccrs.PlateCarree(), ) ax.set_title(panel_title, fontsize=12) shrink_val = 0.90 if view_type == "swath" else 0.78 cb = fig.colorbar(im, ax=ax, orientation="vertical", shrink=shrink_val, pad=0.02) cb.set_label(cbar_label, fontsize=9) fig.suptitle(os.path.basename(datafile) + "\n" + title_suffix, fontsize=12, y=0.98) fig.savefig(output_file, dpi=300, bbox_inches="tight") plt.close(fig) def read_modis_fields(datafile): lat_5km = loadVarFromHDF(datafile, "Latitude") lon_5km = loadVarFromHDF(datafile, "Longitude") CER = clean_fill(loadVarFromHDF(datafile, "Cloud_Effective_Radius") * 0.009999999776482582) CER37 = clean_fill(loadVarFromHDF(datafile, "Cloud_Effective_Radius_37") * 0.009999999776482582) COT = clean_fill(loadVarFromHDF(datafile, "Cloud_Optical_Thickness") * 0.009999999776482582) COT37 = clean_fill(loadVarFromHDF(datafile, "Cloud_Optical_Thickness_37") * 0.009999999776482582) CTH = clean_fill(loadVarFromHDF(datafile, "cloud_top_height_1km") / 1000.0) # MODIS cloud-top temperature is stored with offset + scale in this product. CTT = clean_fill((loadVarFromHDF(datafile, "cloud_top_temperature_1km") + 15000.0) * 0.009999999776482582) SOL_ZEN_5km = loadVarFromHDF(datafile, "Solar_Zenith") * 0.009999999776482582 SEN_ZEN_5km = loadVarFromHDF(datafile, "Sensor_Zenith") * 0.009999999776482582 CF_5km = loadVarFromHDF(datafile, "Cloud_Fraction") * 0.009999999776482582 CM_SPI = loadVarFromHDF(datafile, "Cloud_Mask_SPI") * 0.009999999776482582 CML_flag = loadVarFromHDF(datafile, "Cloud_Multi_Layer_Flag") phaseI_1km = loadVarFromHDF(datafile, "Cloud_Phase_Infrared_1km") CM_SPI1 = CM_SPI[:, :, 0] CM_SPI2 = CM_SPI[:, :, 1] lat_1km = zoom(lat_5km, (CER.shape[0] / lat_5km.shape[0], CER.shape[1] / lat_5km.shape[1])) lon_1km = zoom(lon_5km, (CER.shape[0] / lon_5km.shape[0], CER.shape[1] / lon_5km.shape[1])) # Match the latest uploaded code: repeat 5 km SZA/VZA/CF blocks to 1 km, then crop. ny1, nx1 = CER.shape SOL_ZEN_1km = np.repeat(np.repeat(SOL_ZEN_5km, 5, axis=0), 5, axis=1)[:ny1, :nx1] SEN_ZEN_1km = np.repeat(np.repeat(SEN_ZEN_5km, 5, axis=0), 5, axis=1)[:ny1, :nx1] CF_1km = np.repeat(np.repeat(CF_5km, 5, axis=0), 5, axis=1)[:ny1, :nx1] # Match the latest local-filter code: QA/local filter excludes the final 4 fill columns. lat_1km = trim_last4_columns(lat_1km) lon_1km = trim_last4_columns(lon_1km) CER = trim_last4_columns(CER) CER37 = trim_last4_columns(CER37) COT = trim_last4_columns(COT) COT37 = trim_last4_columns(COT37) CTH = trim_last4_columns(CTH) CTT = trim_last4_columns(CTT) CF_1km = trim_last4_columns(CF_1km) SOL_ZEN_1km = trim_last4_columns(SOL_ZEN_1km) SEN_ZEN_1km = trim_last4_columns(SEN_ZEN_1km) CM_SPI1 = trim_last4_columns(CM_SPI1) CM_SPI2 = trim_last4_columns(CM_SPI2) CML_flag = trim_last4_columns(CML_flag) phaseI_1km = trim_last4_columns(phaseI_1km) # Robust final alignment. Some MODIS arrays can end up as 1350 columns while # QA/local-filter based arrays are 1346 columns after removing the last 4 fill columns. # Crop everything to the smallest common shape before calculating masks/plotting. target_shape = common_2d_shape( lat_1km, lon_1km, CER, CER37, COT, COT37, CTH, CTT, CF_1km, SOL_ZEN_1km, SEN_ZEN_1km, CM_SPI1, CM_SPI2, CML_flag, phaseI_1km, ) lat_1km = crop_to_shape(lat_1km, target_shape) lon_1km = crop_to_shape(lon_1km, target_shape) CER = crop_to_shape(CER, target_shape) CER37 = crop_to_shape(CER37, target_shape) COT = crop_to_shape(COT, target_shape) COT37 = crop_to_shape(COT37, target_shape) CTH = crop_to_shape(CTH, target_shape) CTT = crop_to_shape(CTT, target_shape) CF_1km = crop_to_shape(CF_1km, target_shape) SOL_ZEN_1km = crop_to_shape(SOL_ZEN_1km, target_shape) SEN_ZEN_1km = crop_to_shape(SEN_ZEN_1km, target_shape) CM_SPI1 = crop_to_shape(CM_SPI1, target_shape) CM_SPI2 = crop_to_shape(CM_SPI2, target_shape) CML_flag = crop_to_shape(CML_flag, target_shape) phaseI_1km = crop_to_shape(phaseI_1km, target_shape) LWP = LWP_given_COT_CER(COT, CER) CDNC = CDNC_given_COT_CER(COT, CER) fields = {"LWP": LWP, "CDNC": CDNC, "COT": COT, "CER": CER, "CTH": CTH} filter_inputs = { "CER37": CER37, "COT37": COT37, "SOL_ZEN_1km": SOL_ZEN_1km, "SEN_ZEN_1km": SEN_ZEN_1km, "CM_SPI1": CM_SPI1, "CM_SPI2": CM_SPI2, "CML_flag": CML_flag, "phaseI_1km": phaseI_1km, "CF_1km": CF_1km, "CTT": CTT, } return lat_1km, lon_1km, fields, filter_inputs def apply_latest_filters(datafile, lat_1km, lon_1km, fields, filter_inputs): localFilterMask = computeLocalFilterRemoveMask(datafile, width=local_filter_width) # IMPORTANT: loadVarFromHDF() flips science variables with flipud(). # Flip the QA-based local filter too, as in the latest uploaded code. localFilterMask = np.flipud(localFilterMask) # Force the QA/local-filter mask to exactly match the field/filter array shape. # This avoids 1346-vs-1350 column mismatches from MODIS edge fill columns. ref_shape = fields["CER"].shape localFilterMask = crop_to_shape(localFilterMask, ref_shape) # Base filter from the latest uploaded code, plus the final sampling # constraints you gave in the previous message. latest_mask = ( (filter_inputs["SOL_ZEN_1km"] < SZA_MAX) & (filter_inputs["SEN_ZEN_1km"] < VZA_MAX) & (filter_inputs["CML_flag"] == SINGLE_LAYER_FLAG) & (filter_inputs["phaseI_1km"] == LIQUID_PHASE_IR) & (filter_inputs["CM_SPI1"] < SPI_MAX) & (filter_inputs["CM_SPI2"] < SPI_MAX) & (fields["CER"] >= CER_MIN) & (fields["COT"] >= COT_MIN) & (~localFilterMask) & (fields["COT"] <= COT_MAX) & (filter_inputs["CF_1km"] >= CF_MIN) & (fields["CER"] <= CER_MAX) & (filter_inputs["CTT"] >= CTT_MIN) & (fields["LWP"] >= LWP_MIN) & (fields["LWP"] <= LWP_MAX) & (filter_inputs["CER37"] >= fields["CER"]) ) domain_mask = ( (lat_1km > lat_site - half_width_deg) & (lat_1km < lat_site + half_width_deg) & (lon_1km > lon_site - half_width_deg) & (lon_1km < lon_site + half_width_deg) ) n_domain = int(np.sum(domain_mask)) n_kept_domain = int(np.sum(latest_mask & domain_mask)) print(f" pixels in SGP box kept after updated base filters + final constraints: {n_kept_domain}/{n_domain}") filtered_fields = {} for key, arr in fields.items(): tmp = arr.copy() tmp[~latest_mask] = np.nan filtered_fields[key] = tmp return filtered_fields # ----------------------------- # Main loop # ----------------------------- os.makedirs("MODIS_maps", exist_ok=True) datafiles = sorted(glob("M*.hdf")) print(f"Found {len(datafiles)} MODIS files") for datafile in datafiles: print(datafile) basename = os.path.basename(datafile).replace(".hdf", "") output_unfiltered = f"MODIS_maps/LWP_CDNC_COT_CER_CTH_{basename}_unfiltered_LWPswath_rectPanels_fullData.png" output_filtered = f"MODIS_maps/LWP_CDNC_COT_CER_CTH_{basename}_updatedFilters_finalConstraints_LWPswath_rectPanels_fullData.png" lat_1km, lon_1km, fields, filter_inputs = read_modis_fields(datafile) if not os.path.isfile(output_unfiltered): plot_5panel_figure( datafile, lon_1km, lat_1km, fields, output_unfiltered, title_suffix="Unfiltered pixels, only fill values removed", ) if not os.path.isfile(output_filtered): filtered_fields = apply_latest_filters(datafile, lat_1km, lon_1km, fields, filter_inputs) plot_5panel_figure( datafile, lon_1km, lat_1km, filtered_fields, output_filtered, title_suffix="Updated base filters + COT/CF/CER/CTT/LWP constraints + CER37 ≥ CER", )