From 238921eda39a1af700c0cec7544b66e6e1c70802 Mon Sep 17 00:00:00 2001 From: Jeremy Karst Date: Thu, 16 May 2024 15:43:12 -0400 Subject: [PATCH] All scripts working, initial commit --- ffmpeg_video.py | 59 ++++++ merger.py | 269 ++++++++++++++++++++++++++ merger_FITS.py | 502 ++++++++++++++++++++++++++++++++++++++++++++++++ puller.py | 148 ++++++++++++++ puller_fits.py | 171 +++++++++++++++++ 5 files changed, 1149 insertions(+) create mode 100644 ffmpeg_video.py create mode 100644 merger.py create mode 100644 merger_FITS.py create mode 100644 puller.py create mode 100644 puller_fits.py diff --git a/ffmpeg_video.py b/ffmpeg_video.py new file mode 100644 index 0000000..f88b742 --- /dev/null +++ b/ffmpeg_video.py @@ -0,0 +1,59 @@ +import subprocess +import os +import shutil +import datetime +import time + +import numpy as np +from sortedcontainers import SortedDict +import tqdm + +path_to_images = r"..\composite\goes16" +temp_path = r"..\vid" +output_file = r"..\goes16.mp4" +starttime = time.mktime(datetime.datetime(2024, 1, 1).timetuple()) +stoptime = time.mktime(datetime.datetime(2099, 1, 1).timetuple()) +min_file_size = 200000 # Detect and remove corrupted files by filtering by file size +max_file_size = 500000 + +os.makedirs(temp_path, exist_ok = True) + +files_by_time = SortedDict() +for root, dirs, files in tqdm.tqdm(os.walk(path_to_images), desc="Finding and filtering files."): + for f in files: + if f.endswith('.jpg'): + fpath = os.path.join(root, f) + fsize = os.path.getsize(fpath) + ftime = int(f.replace('-', '.').split('.')[-2]) + if (fsize > min_file_size) and (fsize < max_file_size) and (ftime > starttime) and (ftime < stoptime): + files_by_time[ftime] = fpath + +difftimes = np.diff(files_by_time.keys()) +unique, counts = np.unique(difftimes, return_counts=True) +interval = unique[0] +assert np.sum((unique % interval) > 0) == 0 # Ensure all our timestamps align perfectly with our interval + +prevtime = files_by_time.peekitem(0)[0] - interval +index = 0 +newfpath = "" +for t, f in tqdm.tqdm(files_by_time.items(), desc="Copying files to transcode dir", total = len(files_by_time)): + framejump = (t - prevtime) // interval + if framejump < 60: + for i in range(framejump - 1): # Repeat old frame to fill gaps + shutil.copy(newfpath, f"{os.path.join(temp_path, f"{index+1:06d}.jpg")}") + index += 1 + else: # Skip past intervals that are too large to fill reasonably + print(f"Detected a frame gap of: {framejump}!") + newfpath = f"{os.path.join(temp_path, f"{index+1:06d}.jpg")}" + shutil.copy(f, newfpath) + prevtime = t + index += 1 + +command_line = f'ffmpeg -framerate 60 -pattern_type sequence -i "{os.path.join(temp_path, r"%06d.jpg")}" -c:v libx264 -crf 18 -preset veryfast {output_file}' +print(command_line) + +pipe = subprocess.Popen(command_line, shell=True, stdout=subprocess.PIPE).stdout +output = pipe.read().decode() +pipe.close() + +shutil.rmtree(temp_path) \ No newline at end of file diff --git a/merger.py b/merger.py new file mode 100644 index 0000000..51063cc --- /dev/null +++ b/merger.py @@ -0,0 +1,269 @@ +import os +import time +import datetime +from collections import defaultdict +from multiprocessing import Queue, Process + +import tqdm +from PIL import Image, ImageFont, ImageDraw +import numpy as np +from matplotlib import pyplot as plt + +def bin_ndarray(ndarray, new_shape, operation='mean'): + """ + Bins an ndarray in all axes based on the target shape, by summing or + averaging. + + Number of output dimensions must match number of input dimensions and + new axes must divide old ones. + + Example + ------- + >>> m = np.arange(0,100,1).reshape((10,10)) + >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum') + >>> print(n) + + [[ 22 30 38 46 54] + [102 110 118 126 134] + [182 190 198 206 214] + [262 270 278 286 294] + [342 350 358 366 374]] + + """ + operation = operation.lower() + if not operation in ['sum', 'mean']: + raise ValueError("Operation not supported.") + if ndarray.ndim != len(new_shape): + raise ValueError("Shape mismatch: {} -> {}".format(ndarray.shape, + new_shape)) + compression_pairs = [(d, c//d) for d,c in zip(new_shape, + ndarray.shape)] + flattened = [l for p in compression_pairs for l in p] + ndarray = ndarray.reshape(flattened) + for i in range(len(new_shape)): + op = getattr(ndarray, operation) + ndarray = op(-1*(i+1)) + return ndarray + +def gamma_correct(fun): + def wrapper(*args, **kwargs): + args = list(args) + args[0] = np.power(args[0], 2.2) + args[1] = np.power(args[1], 2.2) + args = tuple(args) + result = fun(*args, **kwargs) + return np.power(result, 1/2.2) + return wrapper + +def clip_color(fun): + def wrapper(*args, **kwargs): + return np.clip(fun(*args, **kwargs), 0.0, 1.0) + return wrapper + +# linear_srgb_matrix = np.array([[0.4124, 0.3576, 0.1805], +# [0.2126, 0.7152, 0.0722], +# [0.0193, 0.1192, 0.9505]]) + +# linear_srgb_matrix_inv = np.array([[ 3.2406, -1.5372, -0.4986], +# [-0.9689, 1.8758, 0.0415], +# [ 0.0557, -0.2040, 1.0570]]) + +# def linear_color_correction(fun): +# def wrapper(*args, **kwargs): +# args = list(args) +# inds = args[0] <= 0.04045 +# ninds = args[0] > 0.04045 +# for i in range(2): +# args[i][inds] = args[i][inds] / 12.92 +# args[i][ninds] = np.power((args[i][ninds] + 0.055) / 1.055, 2.4) +# for x in range(args[i].shape[0]): +# for y in range(args[i].shape[1]): +# args[i][x,y,:] = np.matmul(linear_srgb_matrix, args[i][x,y,:]) +# args = tuple(args) +# result = fun(*args, **kwargs) + +# for x in range(result.shape[0]): +# for y in range(result.shape[1]): +# result[x,y,:] = np.matmul(linear_srgb_matrix_inv, result[x,y,:]) +# inds = result <= 0.0031308 +# ninds = result > 0.0031308 +# result[inds] = result[inds] * 12.92 +# result[ninds] = np.power(result[ninds], 1.0/2.4) * 1.055 - 0.055 +# return result +# return wrapper + + +@gamma_correct +def composite_alpha_over(F, B, alpha_F, alpha_B = 1): + return (F*alpha_F + B*alpha_B*(1-alpha_F)) / (alpha_F + alpha_B*(1-alpha_F)) + +def composite_alpha_blend(F, B, alpha): + return F*alpha + B*(1-alpha) + +def linear_burn(F, B): + burn = F + B - 1 + burn[burn < 0.0] = 0.0 + return burn + +def difference(F, B): + return np.abs(F - B) + +@clip_color +def linear_light(F, B): + result = np.zeros_like(F) + inds = F <= 0.5 + ninds = F > 0.5 + result[inds] = B[inds] + 2.0 * F[inds] - 1 + result[ninds] = 2.0 * (F[ninds] - 0.5) + B[ninds] + return result + +@clip_color +def hard_light(F, B): + result = np.zeros_like(F) + inds = B < 0.5 + ninds = B >= 0.5 + result[inds] = 2 * F[inds] * B[inds] + result[ninds] = 1 - (2*(1 - F[ninds])*(1 - B[ninds])) + return result + +@clip_color +def color_dodge(F, B): + return B / (1.000001 - F) + +@clip_color +def exclusion(F, B): + return F + B - 2*F*B + +@clip_color +def saturation(img, R, G, B): + img[:,:,0] *= R + img[:,:,1] *= G + img[:,:,2] *= B + return img + +@clip_color +def contrast(img, c, b): + return (img - 0.5) * c + 0.5 + b*c + +def generate_composite(work_queue, result_queue): + while True: + try: + job = work_queue.get() + if job is None: + break + files_this_timestamp, timestamp, processed_images_dir = job + filename = f"Composite-{int(timestamp)}.jpg" + filepath = os.path.join(processed_images_dir, filename) + if os.path.isfile(filepath): + result_queue.put(("Exists", timestamp)) + continue + + # image_names = ["094Å", "131Å", "171Å", "195Å", "284Å", "304Å"] + data = [] + for i in range(6): + img = Image.open(files_this_timestamp[i]) + trimmed_img_data = np.array(img)[40:-40,40:-40,:3] # Trim off edges of image to remove text + normalized_img_data = trimmed_img_data / 255.0 # Normalize to float 0.0-1.0 instead of uint8 + srgb_img_data = np.power(normalized_img_data, 2.2) # Gamma correct to sRGB color space + data.append(normalized_img_data) + + # Assemble composite image + composite_image_data = data[4] # Start with (284Å) + # Do a linear burn with 304Å at 95% alpha + composite_image_data = composite_alpha_over(linear_burn(data[5], composite_image_data), composite_image_data, 0.95) + # Do a difference operation with 195Å at 95% alpha + composite_image_data = composite_alpha_over(exclusion(data[3], composite_image_data), composite_image_data, 0.90) + # Do a linear_light layer op with 171Å + composite_image_data = linear_light(data[2], composite_image_data) + # Do a hard_light layer op with 131Å + composite_image_data = composite_alpha_over(hard_light(data[1], composite_image_data), composite_image_data, 0.20) + # Do a color_dodge layer op with 094Å + composite_image_data = composite_alpha_over(color_dodge(data[0], composite_image_data), composite_image_data, 0.25) + # Do an exclusion layer op with 094Å + composite_image_data = composite_alpha_over(exclusion(data[0], composite_image_data), composite_image_data, 0.80) + # Tweak the colors a little + composite_image_data = saturation(composite_image_data, 1.0, 0.95, 1.15) + # Boost contrast + composite_image_data = contrast(composite_image_data, 1.5, 0.15) + + # Now shrink the component images and assemble them alongside the composite. + new_dim = composite_image_data.shape[0] // 3 + # Enlarge the composite to fit the new images + composite_image_data = np.pad(composite_image_data, ((0,0),(new_dim, new_dim),(0,0))) + for i in range(6): + resized = bin_ndarray(data[i], (new_dim, new_dim, 3)) + if i < 3: + composite_image_data[i*new_dim:(i+1)*new_dim, :new_dim, :] = resized + else: + composite_image_data[(i-3)*new_dim:(i-2)*new_dim, -new_dim:, :] = resized + + img = Image.fromarray((255 * composite_image_data).astype('uint8')) + timestring = datetime.datetime.fromtimestamp(timestamp, tz = datetime.UTC).strftime('%Y-%m-%d %H:%M:%S') + ImageDraw.Draw(img).text((655, 15), f"NOAA GOES Sattelite SUVI Composite - {timestring} UTC",(255,255,255), font_size = 24) + img.save(filepath, quality = 90) + result_queue.put(("Created", timestamp)) + + # plt.figure("Composite") + # plt.imshow(composite_image_data) + # plt.show() + # plt.close('all') + except Exception as e: + result_queue.put((e, timestamp)) + + +if __name__ == "__main__": + stored_images_dir = r"..\comp" + processed_images_dir = r"..\composite" + nworkers = 8 + os.makedirs(processed_images_dir, exist_ok=True) + + files_sorted_by_timestamp = defaultdict(list) + for root, dirs, files in os.walk(stored_images_dir): + for f in files: + if f.endswith(".png"): + file_parts = f.split("_") + measurement = file_parts[1] + sattelite = file_parts[2] + measure_end_time = datetime.datetime.strptime(file_parts[4][1:16], "%Y%m%dT%H%M%S") + measure_end_time.replace(tzinfo=datetime.timezone.utc) + measure_end_time = time.mktime(measure_end_time.timetuple()) + files_sorted_by_timestamp[measure_end_time].append(os.path.join(root,f)) + + work_queue = Queue(maxsize = 3) + result_queue = Queue() + workers = [] + for i in range(nworkers): + p = Process(target = generate_composite, args = (work_queue, result_queue), daemon=True) + p.start() + workers.append(p) + + for timestamp in tqdm.tqdm(files_sorted_by_timestamp, desc="Creating Composite Solar Images"): + files_this_timestamp = files_sorted_by_timestamp[timestamp] + files_this_timestamp = sorted(files_this_timestamp) + if not len(files_this_timestamp) == 7: + print(f"Detected a file gap at: {timestamp}") + continue + + work_queue.put((files_this_timestamp, timestamp, processed_images_dir)) + + ncreated = 0 + nexists = 0 + for _ in range(len(files_sorted_by_timestamp)): + result = result_queue.get(5.0) + if result[0] == "Exists": + nexists += 1 + elif result[0] == "Created": + ncreated += 1 + else: + print(f"A worker encountered an exception on job {result[1]}: {result[0]}") + + for _ in range(nworkers): + try: + work_queue.put(None, timeout=1.0) + except: + break + + for w in workers: + w.join(5.0) + + print("Done") \ No newline at end of file diff --git a/merger_FITS.py b/merger_FITS.py new file mode 100644 index 0000000..389005b --- /dev/null +++ b/merger_FITS.py @@ -0,0 +1,502 @@ +import os +import time +import datetime +from collections import defaultdict +from multiprocessing import Queue, Process +import re +import warnings +import queue + +import tqdm +from PIL import Image, ImageDraw +import numpy as np +from matplotlib import pyplot as plt +from astropy.io import fits +import palettable + +def lowpriority(): + """ Set the priority of the process to below-normal.""" + + import sys + try: + sys.getwindowsversion() + except AttributeError: + isWindows = False + else: + isWindows = True + + if isWindows: + # Based on: + # "Recipe 496767: Set Process Priority In Windows" on ActiveState + # http://code.activestate.com/recipes/496767/ + import win32api,win32process,win32con + + pid = win32api.GetCurrentProcessId() + handle = win32api.OpenProcess(win32con.PROCESS_ALL_ACCESS, True, pid) + win32process.SetPriorityClass(handle, win32process.BELOW_NORMAL_PRIORITY_CLASS) + else: + import os + + os.nice(1) + +def bin_ndarray(ndarray, new_shape, operation='mean'): + """ + Bins an ndarray in all axes based on the target shape, by summing or + averaging. + + Number of output dimensions must match number of input dimensions and + new axes must divide old ones. + + Example + ------- + >>> m = np.arange(0,100,1).reshape((10,10)) + >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum') + >>> print(n) + + [[ 22 30 38 46 54] + [102 110 118 126 134] + [182 190 198 206 214] + [262 270 278 286 294] + [342 350 358 366 374]] + + """ + operation = operation.lower() + if not operation in ['sum', 'mean']: + raise ValueError("Operation not supported.") + if ndarray.ndim != len(new_shape): + raise ValueError("Shape mismatch: {} -> {}".format(ndarray.shape, + new_shape)) + compression_pairs = [(d, c//d) for d,c in zip(new_shape, + ndarray.shape)] + flattened = [l for p in compression_pairs for l in p] + ndarray = ndarray.reshape(flattened) + for i in range(len(new_shape)): + op = getattr(ndarray, operation) + ndarray = op(-1*(i+1)) + return ndarray + +def gamma_correct(fun): + def wrapper(*args, **kwargs): + args = list(args) + args[0] = np.power(args[0], 2.2) + args[1] = np.power(args[1], 2.2) + args = tuple(args) + result = fun(*args, **kwargs) + return np.power(result, 1/2.2) + return wrapper + +def clip_color(fun): + def wrapper(*args, **kwargs): + return np.clip(fun(*args, **kwargs), 0.0, 1.0) + return wrapper + +# linear_srgb_matrix = np.array([[0.4124, 0.3576, 0.1805], +# [0.2126, 0.7152, 0.0722], +# [0.0193, 0.1192, 0.9505]]) + +# linear_srgb_matrix_inv = np.array([[ 3.2406, -1.5372, -0.4986], +# [-0.9689, 1.8758, 0.0415], +# [ 0.0557, -0.2040, 1.0570]]) + +# def linear_color_correction(fun): +# def wrapper(*args, **kwargs): +# args = list(args) +# inds = args[0] <= 0.04045 +# ninds = args[0] > 0.04045 +# for i in range(2): +# args[i][inds] = args[i][inds] / 12.92 +# args[i][ninds] = np.power((args[i][ninds] + 0.055) / 1.055, 2.4) +# for x in range(args[i].shape[0]): +# for y in range(args[i].shape[1]): +# args[i][x,y,:] = np.matmul(linear_srgb_matrix, args[i][x,y,:]) +# args = tuple(args) +# result = fun(*args, **kwargs) + +# for x in range(result.shape[0]): +# for y in range(result.shape[1]): +# result[x,y,:] = np.matmul(linear_srgb_matrix_inv, result[x,y,:]) +# inds = result <= 0.0031308 +# ninds = result > 0.0031308 +# result[inds] = result[inds] * 12.92 +# result[ninds] = np.power(result[ninds], 1.0/2.4) * 1.055 - 0.055 +# return result +# return wrapper + + +@gamma_correct +def composite_alpha_over(F, B, alpha_F, alpha_B = 1): + return (F*alpha_F + B*alpha_B*(1-alpha_F)) / (alpha_F + alpha_B*(1-alpha_F)) + +def composite_alpha_blend(F, B, alpha): + return F*alpha + B*(1-alpha) + +def linear_burn(F, B): + burn = F + B - 1 + burn[burn < 0.0] = 0.0 + return burn + +def difference(F, B): + d = np.abs(F - B) + if d.shape[2] == 4: # Preserve alpha of base image + d[:,:,3] = B[:,:,3] + return d + +@clip_color +def linear_light(F, B): + result = np.zeros_like(F) + inds = F <= 0.5 + ninds = F > 0.5 + result[inds] = B[inds] + 2.0 * F[inds] - 1 + result[ninds] = 2.0 * (F[ninds] - 0.5) + B[ninds] + return result + +@clip_color +def hard_light(F, B): + result = np.zeros_like(F) + inds = B < 0.5 + ninds = B >= 0.5 + result[inds] = 2 * F[inds] * B[inds] + result[ninds] = 1 - (2*(1 - F[ninds])*(1 - B[ninds])) + return result + +@clip_color +def color_dodge(F, B): + return B / (1.000001 - F) + +@clip_color +def exclusion(F, B): + d = F + B - 2*F*B + if d.shape[2] == 4: # Preserve alpha of base image + d[:,:,3] = B[:,:,3] + return d + +@clip_color +def saturation(img, R, G, B): + img[:,:,0] *= R + img[:,:,1] *= G + img[:,:,2] *= B + return img + +@clip_color +def contrast(img, c, b): + return (img - 0.5) * c + 0.5 + b*c + +def rgb_to_hsl(img): + r = img[:,:,0] + g = img[:,:,1] + b = img[:,:,2] + cmax = np.copy(r) + cmax[g > cmax] = g[g > cmax] + cmax[b > cmax] = b[b > cmax] + cmin = np.copy(r) + cmin[g < cmin] = g[g < cmin] + cmin[b < cmin] = b[b < cmin] + delta = cmax - cmin + + + # Calc hue + hue = np.zeros_like(r) + inds = cmax == r + with warnings.catch_warnings(): + warnings.filterwarnings('ignore') + hue[inds] = 60 * np.mod((g[inds]-b[inds])/delta[inds], 6) + inds = cmax == g + hue[inds] = 60 * ((b[inds]-r[inds])/delta[inds] + 2) + inds = cmax == b + hue[inds] = 60 * ((r[inds]-g[inds])/delta[inds] + 4) + hue[np.isnan(hue)] = 0 + hue[hue < 0] = hue[hue < 0] + 360 # Make negative hue values positive behind 360 + + # Calc lightness / luminance + luminance = (cmax + cmin) / 2.0 + + # Calc saturation + saturation = np.zeros_like(r) + inds = delta != 0 + saturation[inds] = delta[inds] / (1.0 - np.abs(2 * luminance[inds] - 1.0)) + + # Multiply luminance and saturation by 100 to scale them to the appropriate range (0-100) + saturation *= 100.0 + luminance *= 100.0 + + + if img.shape[2] == 4: # Preserve alpha of original image + return np.stack([hue, saturation, luminance, img[:,:,3]], -1) + else: + return np.stack([hue, saturation, luminance], -1) + +def hsl_to_rgb(img): + hue = img[:,:,0] + saturation = img[:,:,1] + luminance = img[:,:,2] + saturation /= 100.0 + luminance /= 100.0 + c = (1.0 - np.abs(2.0 * luminance - 1.0)) * saturation + x = c * (1.0 - np.abs((hue/60.0) % 2.0 - 1.0)) + m = luminance - c / 2.0 + r = np.zeros_like(hue) + g = np.zeros_like(hue) + b = np.zeros_like(hue) + + inds1 = np.logical_and(0.0 <= hue, hue < 60.0) + r[inds1] = c[inds1] + g[inds1] = x[inds1] + inds2 = np.logical_and(60.0 <= hue, hue < 120.0) + r[inds2] = x[inds2] + g[inds2] = c[inds2] + inds3 = np.logical_and(120.0 <= hue, hue < 180.0) + g[inds3] = c[inds3] + b[inds3] = x[inds3] + inds4 = np.logical_and(180.0 <= hue, hue < 240.0) + g[inds4] = x[inds4] + b[inds4] = c[inds4] + inds5 = np.logical_and(240.0 <= hue, hue < 300.0) + r[inds5] = x[inds5] + b[inds5] = c[inds5] + inds6 = np.logical_and(240.0 <= hue, hue < 300.0) + r[inds6] = c[inds6] + b[inds6] = x[inds6] + + r += m + g += m + b += m + + if img.shape[2] == 4: # Preserve alpha of original image + return np.stack([r, g, b, img[:,:,3]], -1) + else: + return np.stack([r, g, b], -1) + + +def generate_composite(work_queue, result_queue): + image_names = ["094A", "131A", "171A", "195A", "284A", "304A"] + cmaps =[palettable.cmocean.sequential.Ice_5.mpl_colormap, + palettable.cmocean.sequential.Ice_20.mpl_colormap, + palettable.cmocean.sequential.Turbid_5_r.mpl_colormap, + palettable.cmocean.sequential.Turbid_20_r.mpl_colormap, + plt.cm.get_cmap('gist_heat'), + plt.cm.get_cmap('afmhot')] + # Modify colormaps to start at perfect black (when they otherwise start at very dark colors) + for cmi in range(0,4): + for c in ['red','green','blue']: + for i in range(3): + cmaps[cmi]._segmentdata[c][0][i] = 0.0 + # These values are used to map floating point radiance values to colors using the above color maps. + vmins = [0.050, 0.05, 00.100, 00.10, 00.100, 00.1] + vmaxs = [8.000, 8.00, 20.000, 30.00, 40.000, 90.0] + gammas = [0.375, 0.40, 00.425, 00.45, 00.475, 00.5] + while True: + try: + job = work_queue.get() + if job is None: + result_queue.cancel_join_thread() + return + files_this_timestamp, timestamp, processed_images_dir = job + filename = f"Composite-{int(timestamp)}.jpg" + filepath = os.path.join(processed_images_dir, filename) + if os.path.isfile(filepath): + result_queue.put(("Exists", timestamp)) + continue + + base_imgs = [] + for i in range(6): + raw_data = np.flip(fits.getdata(files_this_timestamp[i]), 0) + raw_data[raw_data < 0.0] = 0.0 # Remove non-zero data because it doesn't makes sense (supposed to be std Radiance) + base_imgs.append(cmaps[i](np.clip((raw_data - vmins[i]) / vmaxs[i], 0, 1.0)**gammas[i])) + # plt.figure(image_names[i]) + # plt.imshow(base_imgs[-1]) + + + # plt.figure("Initial Blend") + composite_image_data = composite_alpha_over(base_imgs[4], base_imgs[5], 0.2) + composite_image_data = composite_image_data**0.5 + # plt.imshow(composite_image_data) + + # plt.figure("Linear Burn with 304") + composite_image_data = composite_alpha_over(linear_burn(base_imgs[5], composite_image_data), composite_image_data, 0.60) + # plt.imshow(composite_image_data) + + # plt.figure("Light Ops with mid bands") + mix_img = composite_alpha_over(exclusion(base_imgs[3], composite_image_data), composite_image_data, 0.95) + mix_img = composite_alpha_over(linear_light(base_imgs[2], mix_img), mix_img, 1.0) + mix_img = mix_img**0.75 + composite_image_data = composite_alpha_over(mix_img, composite_image_data, 0.5) + # plt.imshow(composite_image_data) + + # plt.figure("Mix in mid 131") + composite_image_data = composite_alpha_over(base_imgs[2], composite_image_data, 0.25) # Mix in a small amount of 131Å for the nice streamers + # plt.imshow(composite_image_data) + + # plt.figure("Mid diff") + mix_img = difference(base_imgs[3], base_imgs[2]) # 171Å - 131Å + # plt.imshow(mix_img) + + # plt.figure("HSL ops with mid") + comp_hsl = rgb_to_hsl(composite_image_data) + mix_img_hsl = rgb_to_hsl(mix_img) + del mix_img + base_img2_hsl = rgb_to_hsl(base_imgs[2]) + base_img3_hsl = rgb_to_hsl(base_imgs[3]) + comp_hsl = np.copy(comp_hsl) + comp_hsl[:,:,0] += 0.025*mix_img_hsl[:,:,0] # Rotate hue based on mix_1_hsl + del mix_img_hsl + comp_hsl[:,:,0][comp_hsl[:,:,0] > 360.0] -= 360.0 + comp_hsl[:,:,1] -= 0.1*base_img3_hsl[:,:,1] # Reduce saturation based on base_img3 + comp_hsl[:,:,1] = np.clip(comp_hsl[:,:,1], 0.0, 100.0) + comp_hsl[:,:,2] += 0.5*base_img2_hsl[:,:,2] # Boost luminance based on base_img2 + comp_hsl[:,:,1] = np.clip(comp_hsl[:,:,2], 0.0, 100.0) + composite_image_data = hsl_to_rgb(comp_hsl) + del comp_hsl + composite_image_data = saturation(composite_image_data, 1.0, 0.8, 0.1) # Remove some blue and green + # plt.imshow(composite_image_data) + + # plt.figure("Hue shifted 131") + mix_img_hsl = rgb_to_hsl(base_imgs[1]) + mix_img_hsl[:,:,0] -= 50 + mix_img_hsl[:,:,0][mix_img_hsl[:,:,0] < 0.0] += 360.0 + mix_img = hsl_to_rgb(mix_img_hsl) + del mix_img_hsl + # plt.imshow(mix_img) + + # plt.figure("Hard Light with hue shifted 131") + composite_image_data = composite_alpha_over(hard_light(mix_img, composite_image_data), composite_image_data, 0.2) + # plt.imshow(composite_image_data) + + # plt.figure("Adjusted 094") + mix_img = saturation(base_imgs[0], 1.15, 1.2, 1.05) + mix_img = contrast(mix_img, 1.5, 0.0) + # plt.imshow(mix_img) + + # plt.figure("Color Dodge with adjusted 094") + composite_image_data = composite_alpha_over(color_dodge(mix_img, composite_image_data), composite_image_data, 0.5) + # plt.imshow(composite_image_data) + + # plt.figure("Exclusion with adjusted 094") + composite_image_data = composite_alpha_over(exclusion(mix_img, composite_image_data), composite_image_data, 0.65) + # plt.imshow(composite_image_data) + + del mix_img + + # plt.figure("Final Image") + composite_image_data = saturation(composite_image_data, 1.0, 1.1, 1.2) + composite_image_data = contrast(composite_image_data, 1.20, 0.00) + # plt.imshow(composite_image_data) + + # We have our final image + # plt.show() + + # Trim edges of final image so it fits nicely in 1920 + composite_image_data = composite_image_data[64:-64,64:-64,:3] + + # Now shrink the component images and assemble them alongside the composite. + new_dim = composite_image_data.shape[0] // 3 + # Enlarge the composite to fit the new images + composite_image_data = np.pad(composite_image_data, ((0,0),(new_dim, new_dim),(0,0))) + for i in range(6): + img = base_imgs[i][64:-64,64:-64,:3] # Trim edges of image data to fit nicely + img = bin_ndarray(img, (new_dim, new_dim, 3)) # Shrink down to 1/3 for assembly + img = contrast(img, 1.25, 0.0) + xdimoff = i%2 * (composite_image_data.shape[1] - new_dim) + ydimoff = i//2*new_dim + composite_image_data[ydimoff:ydimoff+new_dim, xdimoff:xdimoff+new_dim, :] = img + + img = Image.fromarray((255 * composite_image_data).astype('uint8')) + timestring = datetime.datetime.fromtimestamp(timestamp, tz = datetime.UTC).strftime('%Y-%m-%d %H:%M:%S') + ImageDraw.Draw(img).text((655, 15), f"NOAA GOES Sattelite SUVI Composite - {timestring} UTC",(255,255,255), font_size = 24) + for i in range(6): # Draw component angstrom labels + if i%2 == 0: + xdimtxtoff = 5 + else: + xdimtxtoff = composite_image_data.shape[1] - 44 + ydimtxtoff = i//2*new_dim + new_dim / 2.0 - 8 + ImageDraw.Draw(img).text((xdimtxtoff, ydimtxtoff), image_names[i], font_size = 16) + img.save(filepath, quality = 90) + result_queue.put(("Created", timestamp)) + + except KeyboardInterrupt: + return + except Exception as e: + result_queue.put((e, timestamp)) + + +if __name__ == "__main__": + stored_fits_dirs = [r"..\Data\goes16\l2\data", r"..\Data\goes18\l2\data"] + processed_images_dirs = [r"..\composite\goes16", r"..\composite\goes18"] + + regex_filename = r"dr_suvi-l2-ci\d{3}_g(16|18)_s\S*\.fits" + nworkers = 14 + + # Testing + # stored_fits_dirs = [r"C:\NOAA_SOLAR_IMAGES\fits_test_2024"] + # processed_images_dirs = [r"C:\NOAA_SOLAR_IMAGES\composite"] + + work_queue = Queue(maxsize = nworkers) + result_queue = Queue() + workers = [] + for i in range(nworkers): + p = Process(target = generate_composite, args = (work_queue, result_queue), daemon=True) + p.start() + workers.append(p) + + lowpriority() + ncreated = 0 + nexists = 0 + nfailed = 0 + try: + for stored_fits_dir, processed_images_dir in zip(stored_fits_dirs, processed_images_dirs): + os.makedirs(processed_images_dir, exist_ok=True) + filename_tester = re.compile(regex_filename) + + files_sorted_by_timestamp = defaultdict(list) + found_files = 0 + print(f"Searching for FITS files in: {stored_fits_dir}") + for root, dirs, files in tqdm.tqdm(os.walk(stored_fits_dir), desc="Searching"): + for f in files: + if filename_tester.match(f): + file_parts = f.split("_") + measurement = file_parts[1] + sattelite = file_parts[2] + measure_end_time = datetime.datetime.strptime(file_parts[4][1:16], "%Y%m%dT%H%M%S") + measure_end_time.replace(tzinfo=datetime.timezone.utc) + measure_end_time = time.mktime(measure_end_time.timetuple()) + files_sorted_by_timestamp[measure_end_time].append(os.path.join(root,f)) + found_files += 1 + + print(f"Found {found_files} FITS files. Starting conversion.") + + for timestamp in tqdm.tqdm(files_sorted_by_timestamp, desc="Creating Composite Solar Images"): + # Collect completed jobs and record completion status + while True: + try: + result = result_queue.get_nowait() + if result[0] == "Exists": + nexists += 1 + elif result[0] == "Created": + ncreated += 1 + else: + print(f"A worker encountered an exception on job {result[1]}: {result[0]}") + nfailed += 1 + except queue.Empty: + break + + # Submit new jobs + files_this_timestamp = files_sorted_by_timestamp[timestamp] + files_this_timestamp = sorted(files_this_timestamp) + if not len(files_this_timestamp) == 6: + print(f"Invalid or incomplete sensor records for: {timestamp}") + continue + + work_queue.put((files_this_timestamp, timestamp, processed_images_dir)) + + except KeyboardInterrupt: + print("Finishing current jobs and exiting") + + for _ in range(nworkers): + try: + work_queue.put(None, timeout=10.0) + except: + break + + for w in workers: + w.join(10.0) + + print(f"Created {ncreated} | Already had {nexists} | Failed {nfailed}") \ No newline at end of file diff --git a/puller.py b/puller.py new file mode 100644 index 0000000..e7fe72b --- /dev/null +++ b/puller.py @@ -0,0 +1,148 @@ +import os +import urllib.request +import urllib.parse +import re +import time +import random +import datetime +import json +from threading import Thread +import queue + +import tqdm + +def recursive_find_links(url): + links_regex_pattern = r'(?<=)' # Find href links that do not contain question marks or whitespace + times_regex_pattern = r'(?<=<\/a>)(\s*\d{4}-\d{2}-\d{2} \d{2}:\d{2})(?= )' # Find timestamps in UTC in the format YYYY-MM-DD HH:mm + with urllib.request.urlopen(url) as response: + html_content = response.read().decode('utf-8') + + links = re.findall(links_regex_pattern, html_content) + times = re.findall(times_regex_pattern, html_content) + for i in range(len(times)): + dt = datetime.datetime.strptime(times[i].strip(), "%Y-%m-%d %H:%M") + dt.replace(tzinfo=datetime.timezone.utc) + times[i] = time.mktime(dt.timetuple()) + + if len(links) - 1 == len(times): + links = links[1:] + elif len(links) == len(times): + pass + else: + raise(ValueError) + + for _link, _time in zip(links, times): + if _link.endswith("latest.png"): + continue + elif _link.endswith("/"): + yield from recursive_find_links(urllib.parse.urljoin(url,_link)) + else: + yield urllib.parse.urljoin(url,_link), _time + + +# Fetch image from url and store it to path, retrying on failure +def image_fetch_worker(work_queue, result_queue, attempt_count = 1): + while True: + job = work_queue.get() + if job == None: + return + url, path, t = job + + attempts = 0 + while attempts < attempt_count: + try: + req = urllib.request.Request(url, data=None) + image_data = urllib.request.urlopen(req).read() + os.makedirs(os.path.split(path)[0], exist_ok=True) + open(path, 'wb').write(image_data) + result_queue.put((True, url, t)) + break + except Exception as e: + if hasattr(e, "code") and e.code == 404: # This is expected if the file has been removed from the site (at least for swpc.noaa.gov) + attempts += attempt_count + elif attempts == 0: + print(f"\nA problem occurred on image: {url} | {e}") + time.sleep(1 + random.random()) + attempts += 1 + if (attempt_count > 1) and (attempts == attempt_count): + print(f'\nAfter {attempt_count} retries, could not fetch: {url}') + result_queue.put((False, url, t)) + break + + +if __name__ == "__main__": + directory_url = r"https://services.swpc.noaa.gov/images/animations/suvi/" + stored_images_dir = r"..\suvi" + file_database_path = r"..\file_database.json" + fetch_interval = 30*60 + nworkers = 8 + + file_info_cache = {} + try: + print(f"Attempting to load file records from cache: {file_database_path}") + with open(file_database_path, 'r') as f: + file_info_cache = json.loads(f.read()) + print(f"File records loaded from cache: {len(file_info_cache)} records found.") + except Exception as e: + print(f"Load failed, starting with empty cache") + file_info_cache = {} + + work_queue = queue.Queue(maxsize=nworkers) + result_queue = queue.Queue() + workers = [] + for _ in range(nworkers): + t = Thread(target=image_fetch_worker, args=(work_queue, result_queue), daemon=True) + t.start() + workers.append(t) + try: + while True: + fetched_image_count = 0 + already_had_image_count = 0 + failed_image_count = 0 + urllen = len(directory_url) + for l, t in tqdm.tqdm(recursive_find_links(directory_url), desc="Downloading files"): + # Collect complete work and record it + while True: + try: + r_success, r_url, r_t = result_queue.get_nowait() + if r_success: + fetched_image_count += 1 + file_info_cache[r_url] = r_t + else: + failed_image_count += 1 + except queue.Empty: + break + # If we dont have the file or the file at the link is newer than the one we previously fetched + if (not (l in file_info_cache)) or t > file_info_cache[l]: + file_portion_of_link = l[urllen:] + filepath = os.path.join(stored_images_dir, file_portion_of_link.replace("/", os.sep)) + if os.path.exists(filepath): + if l in file_info_cache: # If we have record of this file, it must be out of date, rename it and download the new version. + os.rename(filepath, filepath + f"_{int(file_info_cache[l])}") + else: # If we have no record of this file, update the file info cache and don't redownload + file_info_cache[l] = t + already_had_image_count += 1 + continue + work_queue.put((l, filepath, t)) + else: + already_had_image_count += 1 + with open(file_database_path, 'w') as f: + f.write(json.dumps(file_info_cache)) + print(f"Downloaded {fetched_image_count} | Already had {already_had_image_count} | Failed {failed_image_count}") + time.sleep(fetch_interval) + except KeyboardInterrupt: + print("Saving file database and shutting down.") + + for _ in range(nworkers): + try: + work_queue.put(None, timeout=1.0) + except: + break + + for w in workers: + w.join(5.0) + + with open(file_database_path, 'w') as f: + f.write(json.dumps(file_info_cache)) + + print(f"Downloaded {fetched_image_count} | Already had {already_had_image_count} | Failed {failed_image_count}") \ No newline at end of file diff --git a/puller_fits.py b/puller_fits.py new file mode 100644 index 0000000..24b2d0c --- /dev/null +++ b/puller_fits.py @@ -0,0 +1,171 @@ +import os +import urllib.request +import urllib.parse +import re +import time +import random +import datetime +import json +from threading import Thread +import queue + +import tqdm + +directory_url = r"https://data.ngdc.noaa.gov/platforms/solar-space-observing-satellites/goes/" +stored_images_dir = r"..\Data" +# directory_url = r"https://data.ngdc.noaa.gov/platforms/solar-space-observing-satellites/goes/goes16/l2/data/suvi-l2-ci094/" +# stored_images_dir = r"Z:\NOAA GOES Data\Data\goes16\l2\suvi-l2-ci094" +ignore_folder_names = ["l1b", "goes17", "2017", "2018", "2019", "2020", "2021", "2022"] +file_database_path = r"..\file_database.json" +fetch_interval = 60*60 +nworkers = 16 +randomize_order = True + +def recursive_find_links(url, attempt_count = 3, randomize_order = True): + links_regex_pattern = r'(?<=.*\d{4}-\d{2}-\d{2} \d{2}:\d{2})' # Find href links that do not contain question marks or whitespace + times_regex_pattern = r'(?<=)(\d{4}-\d{2}-\d{2} \d{2}:\d{2})(?= )' # Find timestamps in UTC in the format YYYY-MM-DD HH:mm + with urllib.request.urlopen(url) as response: + attempts = 0 + html_content = None + while attempts < attempt_count: + try: + html_content = response.read().decode('utf-8') + break + except Exception as e: + print(f"Exception while fetching links: {e}") + time.sleep(1 + random.random()) + attempts += 1 + if (attempt_count > 1) and (attempts == attempt_count): + print(f'\nAfter {attempt_count} retries, could not fetch: {url}') + return + + links = re.findall(links_regex_pattern, html_content) + times = re.findall(times_regex_pattern, html_content) + for i in range(len(times)): + dt = datetime.datetime.strptime(times[i].strip(), "%Y-%m-%d %H:%M") + dt.replace(tzinfo=datetime.timezone.utc) + times[i] = time.mktime(dt.timetuple()) + + if len(links) - 1 == len(times): + links = links[1:] + elif len(links) == len(times): + pass + else: + raise(ValueError) + + results = list(zip(links, times)) + if randomize_order: + random.shuffle(results) + for _link, _time in results: + if _link.endswith("/"): + if _link.split(r"/")[-2] in ignore_folder_names: + continue + else: + yield from recursive_find_links(urllib.parse.urljoin(url,_link), randomize_order) + else: + yield urllib.parse.urljoin(url,_link), _time + + +# Fetch image from url and store it to path, retrying on failure +def image_fetch_worker(work_queue, result_queue, attempt_count = 1): + while True: + job = work_queue.get() + if job == None: + return + url, path, t = job + + attempts = 0 + while attempts < attempt_count: + try: + req = urllib.request.Request(url, data=None) + image_data = urllib.request.urlopen(req).read() + os.makedirs(os.path.split(path)[0], exist_ok=True) + open(path, 'wb').write(image_data) + result_queue.put((True, url, t)) + break + except Exception as e: + if hasattr(e, "code") and e.code == 404: # This is expected if the file has been removed from the site (at least for swpc.noaa.gov) + attempts += attempt_count + elif attempts == 0: + print(f"\nA problem occurred on image: {url} | {e}") + time.sleep(1 + random.random()) + attempts += 1 + if (attempt_count > 1) and (attempts == attempt_count): + print(f'\nAfter {attempt_count} retries, could not fetch: {url}') + result_queue.put((False, url, t)) + break + + +if __name__ == "__main__": + file_info_cache = {} + try: + print(f"Attempting to load file records from cache: {file_database_path}") + with open(file_database_path, 'r') as f: + file_info_cache = json.loads(f.read()) + print(f"File records loaded from cache: {len(file_info_cache)} records found.") + except Exception as e: + print(f"Load failed, starting with empty cache") + file_info_cache = {} + + work_queue = queue.Queue(maxsize=nworkers) + result_queue = queue.Queue() + workers = [] + for _ in range(nworkers): + t = Thread(target=image_fetch_worker, args=(work_queue, result_queue), daemon=True) + t.start() + workers.append(t) + try: + while True: + fetched_image_count = 0 + already_had_image_count = 0 + failed_image_count = 0 + urllen = len(directory_url) + for l, t in tqdm.tqdm(recursive_find_links(directory_url, randomize_order), desc="Downloading files"): + # Collect completed jobs and record completion status + while True: + try: + r_success, r_url, r_t = result_queue.get_nowait() + if r_success: + fetched_image_count += 1 + file_info_cache[r_url] = r_t + else: + failed_image_count += 1 + except queue.Empty: + break + # If we dont have the file or the file at the link is newer than the one we previously fetched + if (not (l in file_info_cache)) or t > file_info_cache[l]: + file_portion_of_link = l[urllen:] + filepath = os.path.join(stored_images_dir, file_portion_of_link.replace("/", os.sep)) + if os.path.exists(filepath): + if l in file_info_cache: # If we have record of this file, it must be out of date, rename it and download the new version. + new_name = filepath + f"_{int(file_info_cache[l])}" + if os.path.exists(new_name): + os.remove(new_name) + os.rename(filepath, new_name) + else: # If we have no record of this file, update the file info cache and don't redownload + file_info_cache[l] = t + already_had_image_count += 1 + continue + work_queue.put((l, filepath, t)) + else: + already_had_image_count += 1 + with open(file_database_path, 'w') as f: + f.write(json.dumps(file_info_cache)) + print(f"Downloaded {fetched_image_count} | Already had {already_had_image_count} | Failed {failed_image_count}") + time.sleep(fetch_interval) + except KeyboardInterrupt: + print("Saving file database and shutting down.") + + for _ in range(nworkers): + try: + work_queue.put(None, timeout=5.0) + except: + break + + for w in workers: + w.join(5.0) + + with open(file_database_path, 'w') as f: + f.write(json.dumps(file_info_cache)) + + print(f"Downloaded {fetched_image_count} | Already had {already_had_image_count} | Failed {failed_image_count}") \ No newline at end of file