From 5a30bef3937ee856c0f56b7abf97e947296296a2 Mon Sep 17 00:00:00 2001 From: Hayley Clevenger Date: Thu, 25 Jul 2024 19:57:35 -0400 Subject: [PATCH 1/3] adding some new plotting scripts --- hc_processing_scripts/alt_cut_plotting.py | 65 +++++++++++++++++++++++ hc_processing_scripts/lat_cut_plotting.py | 65 +++++++++++++++++++++++ hc_processing_scripts/lon_cut_plotting.py | 65 +++++++++++++++++++++++ 3 files changed, 195 insertions(+) create mode 100644 hc_processing_scripts/alt_cut_plotting.py create mode 100644 hc_processing_scripts/lat_cut_plotting.py create mode 100644 hc_processing_scripts/lon_cut_plotting.py diff --git a/hc_processing_scripts/alt_cut_plotting.py b/hc_processing_scripts/alt_cut_plotting.py new file mode 100644 index 0000000..a421ffa --- /dev/null +++ b/hc_processing_scripts/alt_cut_plotting.py @@ -0,0 +1,65 @@ +import gemini3d.read as read +from gemini3d.grid.gridmodeldata import model2geogcoords +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.cm import ScalarMappable +from matplotlib.colors import Normalize + +# File setup +direc = '/Users/clevenger/Simulations/poker_asi_mar19/' # path to directory containing simulation output files +outdir = '/Users/clevenger/Simulations/poker_asi_mar19/postprocessing/' # path to save plot outputs + +# File reading stuff +cfg = read.config(direc) +xg = read.grid(direc) +parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 +dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index + +lalt = 128; llon = 128; llat = 128; # resolutions to plot model outputs (higher=better) + +print("Sampling in geographic coords...") +galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) + +print("Creating altitude slices...") +altitudes = [110e3, 310e3, 510e3, 710e3] +alt_indices = [np.argmin(abs(galti - alt)) for alt in altitudes] + +# Lock latitude and longitude into place +ilon = np.argmin(abs(gloni - 212)) +ilat = np.argmin(abs(glati - 65)) + +# Creating 3D figure space to put 2D maps onto +fig = plt.figure(figsize=(15, 8)) +ax = fig.add_subplot(111, projection='3d') +ax.grid(False) + +# Normalize plot color/quantities so all of the 2D maps can share the same colorbar +norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) +mappable = ScalarMappable(norm=norm, cmap='viridis') + +print("Making altitude cuts...") +for alt_idx in alt_indices: + X, Y = np.meshgrid(gloni, glati) # latitude and longitude + Z = np.full_like(X, galti[alt_idx]) # altitude + C = parmgi[alt_idx, :, :].T # extract data at fixed altitude + C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping + + print("Plotting each altitude cut...") + ax.plot_surface(X, Y, Z, facecolors=plt.cm.viridis(C_normalized), rstride=1, cstride=1, edgecolor='none') + +print("Making combined altitude slice figure...") +# Add colorbar + labels +cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) +cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)' +ax.set_xlabel('Latitude (deg)') +ax.set_ylabel('Longitude (deg)') +ax.set_zlabel('Altitude (m)') +ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want +ax.set_title('Electron Density at various altitudes') # change for whatever parm + +print("Saving figure...") +# Save + show plots +filename = 'ne_alt_cuts.png' # adjust filename for plot being saved +plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution + +plt.show() \ No newline at end of file diff --git a/hc_processing_scripts/lat_cut_plotting.py b/hc_processing_scripts/lat_cut_plotting.py new file mode 100644 index 0000000..08fc688 --- /dev/null +++ b/hc_processing_scripts/lat_cut_plotting.py @@ -0,0 +1,65 @@ +import gemini3d.read as read +from gemini3d.grid.gridmodeldata import model2geogcoords +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.cm import ScalarMappable +from matplotlib.colors import Normalize + +# File setup +direc = '/Users/clevenger/Simulations/poker_asi_mar19/' # path to directory containing simulation output files +outdir = '/Users/clevenger/Simulations/poker_asi_mar19/postprocessing/' # path to save plot outputs + +# File reading stuff +cfg = read.config(direc) +xg = read.grid(direc) +parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 +dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index + +lalt = 128; llon = 128; llat = 128; # resolutions to plot model outputs (higher=better) + +print("Sampling in geographic coords...") +galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) + +print("Creating latitude slices...") +latitudes = [61, 63, 65, 67] # add/change latitudes as needed +lat_indices = [np.argmin(abs(glati - lat)) for lat in latitudes] + +# Lock longitude and altitude into place +ilon = np.argmin(abs(gloni - 212)) +ialt = np.argmin(abs(galti - 110)) + +# Creating 3D figure space to put 2D maps onto +fig = plt.figure(figsize=(10, 8)) +ax = fig.add_subplot(111, projection='3d') +ax.grid(False) + +# Normalize plot color/quantities so all of the 2D maps can share the same colorbar +norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) +mappable = ScalarMappable(norm=norm, cmap='viridis') + +print("Making latitude cuts...") +for lat_idx in lat_indices: + X, Z = np.meshgrid(gloni, galti) # longitude and altitude + Y = np.full_like(X, glati[lat_idx]) # latitude + C = parmgi[:, :, lat_idx].T # extract data at fixed latitude + C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping + + print("Plotting each latitude cut...") + ax.plot_surface(X, Y, Z, facecolors=plt.cm.viridis(C_normalized), rstride=1, cstride=1, edgecolor='none') + +print("Making combined altitude slice figure...") +# Add colorbar + labels +cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) +cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)' +ax.set_xlabel('Longitude (deg)') +ax.set_ylabel('Latitude (deg)') +ax.set_zlabel('Altitude (m)') +ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want +ax.set_title('Electron Density at various latitudes') # change for whatever parm + +print("Saving figure...") +# Save + show plots +filename = 'ne_lat_cuts.png' # adjust filename for plot being saved +plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution + +plt.show() diff --git a/hc_processing_scripts/lon_cut_plotting.py b/hc_processing_scripts/lon_cut_plotting.py new file mode 100644 index 0000000..9b67496 --- /dev/null +++ b/hc_processing_scripts/lon_cut_plotting.py @@ -0,0 +1,65 @@ +import gemini3d.read as read +from gemini3d.grid.gridmodeldata import model2geogcoords +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.cm import ScalarMappable +from matplotlib.colors import Normalize + +# File setup +direc = '/Users/clevenger/Simulations/poker_asi_mar19/' # path to directory containing simulation output files +outdir = '/Users/clevenger/Simulations/poker_asi_mar19/postprocessing/' # path to save plot outputs + +# File reading stuff +cfg = read.config(direc) +xg = read.grid(direc) +parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 +dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index + +lalt = 128; llon = 128; llat = 128; # resolutions to plot model outputs (higher=better) + +print("Sampling in geographic coords...") +galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) + +print("Creating longitude slices...") +longitudes = [200, 205, 210, 215, 220] # add/change latitudes as needed +lon_indices = [np.argmin(abs(gloni - lon)) for lon in longitudes] + +# Lock latitude and altitude into place +ilat = np.argmin(abs(glati - 65)) +ialt = np.argmin(abs(galti - 110)) + +# Creating 3D figure space to put 2D maps onto +fig = plt.figure(figsize=(10, 8)) +ax = fig.add_subplot(111, projection='3d') +ax.grid(False) + +# Normalize plot color/quantities so all of the 2D maps can share the same colorbar +norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) +mappable = ScalarMappable(norm=norm, cmap='viridis') + +print("Making longitude cuts...") +for lon_idx in lon_indices: + X, Z = np.meshgrid(glati, galti) # latitude and altitude + Y = np.full_like(X, gloni[lon_idx]) # longitude + C = parmgi[:, lon_idx, :].T # extract data at fixed longitude + C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping + + print("Plotting each longitude cut...") + ax.plot_surface(X, Y, Z, facecolors=plt.cm.viridis(C_normalized), rstride=1, cstride=1, edgecolor='none') + +print("Making combined longitude slice figure...") +# Add colorbar + labels +cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) +cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)' +ax.set_xlabel('Latitude (deg)') +ax.set_ylabel('Longitude (deg)') +ax.set_zlabel('Altitude (m)') +ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want +ax.set_title('Electron Density at various latitudes') # change for whatever parm + +print("Saving figure...") +# Save + show plots +filename = 'ne_lon_cuts.png' # adjust filename for plot being saved +plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution + +plt.show() \ No newline at end of file From ed669fe5ab8066df6ca058dc95a3bcba8793633d Mon Sep 17 00:00:00 2001 From: Hayley Clevenger Date: Fri, 26 Dec 2025 21:29:33 -0500 Subject: [PATCH 2/3] adding particles changes for Q inputs from dasc inversions (and a few plotting scripts) --- ..._postprocessing_plasma_params_automated.py | 76 +++++++++++++++++++ .../alt_cut_plotting.py | 16 ++-- .../lat_cut_plotting.py | 0 .../lon_cut_plotting.py | 0 src/gemini3d/particles/__init__.py | 18 ++++- src/gemini3d/particles/gaussian2d.py | 4 +- 6 files changed, 101 insertions(+), 13 deletions(-) create mode 100644 hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py rename {hc_processing_scripts => hz_processing_scripts}/alt_cut_plotting.py (82%) rename {hc_processing_scripts => hz_processing_scripts}/lat_cut_plotting.py (100%) rename {hc_processing_scripts => hz_processing_scripts}/lon_cut_plotting.py (100%) diff --git a/hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py b/hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py new file mode 100644 index 0000000..e3200ce --- /dev/null +++ b/hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py @@ -0,0 +1,76 @@ +import gemini3d.read as read +from gemini3d.grid.gridmodeldata import model2geogcoords +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.cm import ScalarMappable +from matplotlib.colors import Normalize + +# File setup +direc = '/Users/clevenger/Simulations/agu24/fac_precip/' # path to directory containing simulation output files +outdir = '/Users/clevenger/Simulations/agu24/fac_precip/postprocessing/' # path to save plot outputs + +# File reading stuff +cfg = read.config(direc) +xg = read.grid(direc) +parm = "v1" # choose from ne, Ti, Te, v1, v2, v3, J1, J2, or J3 + +lalt = 256 +llon = 256 +llat = 256 # resolutions to plot model outputs (higher = better) + +# User-specified time range +start_index = 60 # Adjust as needed +end_index = 120 # Adjust as needed + +# Loop over the desired range of time indices +for time_idx in range(start_index, end_index + 1): + time_value = cfg["time"][time_idx] # Get the current time value + print(f"Processing time index {time_idx} with time {time_value}...") + + # Read the data for the current time + dat = read.frame(direc, time_value, var=parm) + galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) + + # Altitude cuts + altitudes = [110e3, 310e3, 510e3, 710e3] + alt_indices = [np.argmin(abs(galti - alt)) for alt in altitudes] + + # Lock latitude and longitude into place + ilon = np.argmin(abs(gloni - 212)) + ilat = np.argmin(abs(glati - 65)) + + # Create 3D figure + fig = plt.figure(figsize=(15, 8)) + ax = fig.add_subplot(111, projection='3d') + ax.grid(False) + + # Normalize plot color/quantities + norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) + mappable = ScalarMappable(norm=norm, cmap='bwr') + + # Plot altitude cuts + for alt_idx in alt_indices: + X, Y = np.meshgrid(gloni, glati) # latitude and longitude + Z = np.full_like(X, galti[alt_idx]) # altitude + C = parmgi[alt_idx, :, :].T # extract data at fixed altitude + C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping + + print("Plotting each altitude cut...") + ax.plot_surface(X, Y, Z, facecolors=plt.cm.bwr(C_normalized), rstride=1, cstride=1, edgecolor='none') + + print("Making combined altitude slice figure...") + # Add colorbar + labels + cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) + cbar.set_label('v$_{1} (m/s)$') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)', 'ne (m$^{-3}$)' + ax.set_xlabel('Latitude (deg)') + ax.set_ylabel('Longitude (deg)') + ax.set_zlabel('Altitude (m)') + ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want + ax.set_title(f'Line of Sight Velocity at Various Altitudes\nTime: {time_value.strftime("%Y-%m-%d %H:%M:%S")}') # change for whatever parm + + # Save the figure + filename = f'v1_alt_cut_{time_value.strftime("%Y%m%d_%H%M%S")}.png' # Format time for filename + plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') + print(f"Saved figure: {outdir}/{filename}") + + plt.close(fig) # Close the figure to free up memory diff --git a/hc_processing_scripts/alt_cut_plotting.py b/hz_processing_scripts/alt_cut_plotting.py similarity index 82% rename from hc_processing_scripts/alt_cut_plotting.py rename to hz_processing_scripts/alt_cut_plotting.py index a421ffa..d236f73 100644 --- a/hc_processing_scripts/alt_cut_plotting.py +++ b/hz_processing_scripts/alt_cut_plotting.py @@ -6,16 +6,16 @@ from matplotlib.colors import Normalize # File setup -direc = '/Users/clevenger/Simulations/poker_asi_mar19/' # path to directory containing simulation output files -outdir = '/Users/clevenger/Simulations/poker_asi_mar19/postprocessing/' # path to save plot outputs +direc = '/Users/clevenger/Simulations/paper01_event01/' # path to directory containing simulation output files +outdir = '/Users/clevenger/Projects/paper01/events/20230227/gemini_postprocessing/' # path to save plot outputs # File reading stuff cfg = read.config(direc) xg = read.grid(direc) parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 -dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index +dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index START AT 61 -lalt = 128; llon = 128; llat = 128; # resolutions to plot model outputs (higher=better) +lalt = 256; llon = 256; llat = 256; # resolutions to plot model outputs (higher=better) print("Sampling in geographic coords...") galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) @@ -35,7 +35,7 @@ # Normalize plot color/quantities so all of the 2D maps can share the same colorbar norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) -mappable = ScalarMappable(norm=norm, cmap='viridis') +mappable = ScalarMappable(norm=norm, cmap='gnuplot2') print("Making altitude cuts...") for alt_idx in alt_indices: @@ -45,12 +45,12 @@ C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping print("Plotting each altitude cut...") - ax.plot_surface(X, Y, Z, facecolors=plt.cm.viridis(C_normalized), rstride=1, cstride=1, edgecolor='none') + ax.plot_surface(X, Y, Z, facecolors=plt.cm.gnuplot2(C_normalized), rstride=1, cstride=1, edgecolor='none') print("Making combined altitude slice figure...") # Add colorbar + labels cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) -cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)' +cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)', 'ne (m$^{-3}$)' ax.set_xlabel('Latitude (deg)') ax.set_ylabel('Longitude (deg)') ax.set_zlabel('Altitude (m)') @@ -62,4 +62,4 @@ filename = 'ne_alt_cuts.png' # adjust filename for plot being saved plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution -plt.show() \ No newline at end of file +plt.show() diff --git a/hc_processing_scripts/lat_cut_plotting.py b/hz_processing_scripts/lat_cut_plotting.py similarity index 100% rename from hc_processing_scripts/lat_cut_plotting.py rename to hz_processing_scripts/lat_cut_plotting.py diff --git a/hc_processing_scripts/lon_cut_plotting.py b/hz_processing_scripts/lon_cut_plotting.py similarity index 100% rename from hc_processing_scripts/lon_cut_plotting.py rename to hz_processing_scripts/lon_cut_plotting.py diff --git a/src/gemini3d/particles/__init__.py b/src/gemini3d/particles/__init__.py index 40de7c3..da073dc 100644 --- a/src/gemini3d/particles/__init__.py +++ b/src/gemini3d/particles/__init__.py @@ -7,6 +7,8 @@ from .core import get_times from .grid import precip_grid +# OUTPUT ARGUMENT NEEDS TO BE MADE OPTIONAL TO RUN 'OLD' EXAMPLES!!!!! + # this is loaded dynamically via str2func from .gaussian2d import gaussian2d @@ -52,17 +54,25 @@ def particles_BCs(cfg: dict[str, T.Any], xg: dict[str, T.Any]): else: Qfunc = str2func("gemini3d.particles.gaussian2d") - Qtmp = Qfunc(pg, cfg["Qprecip"], cfg["Qprecip_background"]) - + Qtmp, E0_temp = Qfunc(pg, cfg["Qprecip"], cfg["Qprecip_background"]) + print("i_on: ", i_on) + print("i_off: ", i_off) + print("Qtmp shape: ", Qtmp.shape) + print("E0tmp shape: ", E0_temp.shape) + for i in range(i_on, i_off): if Qtmp.ndim == 3: - pg["Q"][i, :, :] = Qtmp[i, :, :] + pg["Q"][i, :, :] = Qtmp[i, :, :] # time, lon, lat else: pg["Q"][i, :, :] = Qtmp - pg["E0"][i, :, :] = cfg["E0precip"] + pg["E0"][i, :, :] = E0_temp assert np.isfinite(pg["Q"]).all(), "Q flux must be finite" assert (pg["Q"] >= 0).all(), "Q flux must be non-negative" + if E0 > 0 + assert (pg["E0"] >=0).all(), "E0 characteristic energy must be non-negative" + else + continue # %% CONVERT THE ENERGY TO EV # E0 = max(E0,0.100); diff --git a/src/gemini3d/particles/gaussian2d.py b/src/gemini3d/particles/gaussian2d.py index e48b3a2..3ef1b60 100644 --- a/src/gemini3d/particles/gaussian2d.py +++ b/src/gemini3d/particles/gaussian2d.py @@ -31,5 +31,7 @@ def gaussian2d(pg, Qpeak: float, Qbackground: float): raise LookupError("precipation must be defined in latitude, longitude or both") Q[Q < Qbackground] = Qbackground + + E0 = cfg["E0precip"] - return Q + return Q, E0 From a14c1503b39d66e4a735589203001a9632a9087f Mon Sep 17 00:00:00 2001 From: Hayley Clevenger Date: Fri, 26 Dec 2025 21:40:10 -0500 Subject: [PATCH 3/3] taking out plotting scripts to directly link to issue #23 without add-ons --- ..._postprocessing_plasma_params_automated.py | 76 ------------------- hz_processing_scripts/alt_cut_plotting.py | 65 ---------------- hz_processing_scripts/lat_cut_plotting.py | 65 ---------------- hz_processing_scripts/lon_cut_plotting.py | 65 ---------------- 4 files changed, 271 deletions(-) delete mode 100644 hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py delete mode 100644 hz_processing_scripts/alt_cut_plotting.py delete mode 100644 hz_processing_scripts/lat_cut_plotting.py delete mode 100644 hz_processing_scripts/lon_cut_plotting.py diff --git a/hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py b/hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py deleted file mode 100644 index e3200ce..0000000 --- a/hz_processing_scripts/WIP/gemini_postprocessing_plasma_params_automated.py +++ /dev/null @@ -1,76 +0,0 @@ -import gemini3d.read as read -from gemini3d.grid.gridmodeldata import model2geogcoords -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.cm import ScalarMappable -from matplotlib.colors import Normalize - -# File setup -direc = '/Users/clevenger/Simulations/agu24/fac_precip/' # path to directory containing simulation output files -outdir = '/Users/clevenger/Simulations/agu24/fac_precip/postprocessing/' # path to save plot outputs - -# File reading stuff -cfg = read.config(direc) -xg = read.grid(direc) -parm = "v1" # choose from ne, Ti, Te, v1, v2, v3, J1, J2, or J3 - -lalt = 256 -llon = 256 -llat = 256 # resolutions to plot model outputs (higher = better) - -# User-specified time range -start_index = 60 # Adjust as needed -end_index = 120 # Adjust as needed - -# Loop over the desired range of time indices -for time_idx in range(start_index, end_index + 1): - time_value = cfg["time"][time_idx] # Get the current time value - print(f"Processing time index {time_idx} with time {time_value}...") - - # Read the data for the current time - dat = read.frame(direc, time_value, var=parm) - galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) - - # Altitude cuts - altitudes = [110e3, 310e3, 510e3, 710e3] - alt_indices = [np.argmin(abs(galti - alt)) for alt in altitudes] - - # Lock latitude and longitude into place - ilon = np.argmin(abs(gloni - 212)) - ilat = np.argmin(abs(glati - 65)) - - # Create 3D figure - fig = plt.figure(figsize=(15, 8)) - ax = fig.add_subplot(111, projection='3d') - ax.grid(False) - - # Normalize plot color/quantities - norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) - mappable = ScalarMappable(norm=norm, cmap='bwr') - - # Plot altitude cuts - for alt_idx in alt_indices: - X, Y = np.meshgrid(gloni, glati) # latitude and longitude - Z = np.full_like(X, galti[alt_idx]) # altitude - C = parmgi[alt_idx, :, :].T # extract data at fixed altitude - C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping - - print("Plotting each altitude cut...") - ax.plot_surface(X, Y, Z, facecolors=plt.cm.bwr(C_normalized), rstride=1, cstride=1, edgecolor='none') - - print("Making combined altitude slice figure...") - # Add colorbar + labels - cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) - cbar.set_label('v$_{1} (m/s)$') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)', 'ne (m$^{-3}$)' - ax.set_xlabel('Latitude (deg)') - ax.set_ylabel('Longitude (deg)') - ax.set_zlabel('Altitude (m)') - ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want - ax.set_title(f'Line of Sight Velocity at Various Altitudes\nTime: {time_value.strftime("%Y-%m-%d %H:%M:%S")}') # change for whatever parm - - # Save the figure - filename = f'v1_alt_cut_{time_value.strftime("%Y%m%d_%H%M%S")}.png' # Format time for filename - plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') - print(f"Saved figure: {outdir}/{filename}") - - plt.close(fig) # Close the figure to free up memory diff --git a/hz_processing_scripts/alt_cut_plotting.py b/hz_processing_scripts/alt_cut_plotting.py deleted file mode 100644 index d236f73..0000000 --- a/hz_processing_scripts/alt_cut_plotting.py +++ /dev/null @@ -1,65 +0,0 @@ -import gemini3d.read as read -from gemini3d.grid.gridmodeldata import model2geogcoords -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.cm import ScalarMappable -from matplotlib.colors import Normalize - -# File setup -direc = '/Users/clevenger/Simulations/paper01_event01/' # path to directory containing simulation output files -outdir = '/Users/clevenger/Projects/paper01/events/20230227/gemini_postprocessing/' # path to save plot outputs - -# File reading stuff -cfg = read.config(direc) -xg = read.grid(direc) -parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 -dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index START AT 61 - -lalt = 256; llon = 256; llat = 256; # resolutions to plot model outputs (higher=better) - -print("Sampling in geographic coords...") -galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) - -print("Creating altitude slices...") -altitudes = [110e3, 310e3, 510e3, 710e3] -alt_indices = [np.argmin(abs(galti - alt)) for alt in altitudes] - -# Lock latitude and longitude into place -ilon = np.argmin(abs(gloni - 212)) -ilat = np.argmin(abs(glati - 65)) - -# Creating 3D figure space to put 2D maps onto -fig = plt.figure(figsize=(15, 8)) -ax = fig.add_subplot(111, projection='3d') -ax.grid(False) - -# Normalize plot color/quantities so all of the 2D maps can share the same colorbar -norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) -mappable = ScalarMappable(norm=norm, cmap='gnuplot2') - -print("Making altitude cuts...") -for alt_idx in alt_indices: - X, Y = np.meshgrid(gloni, glati) # latitude and longitude - Z = np.full_like(X, galti[alt_idx]) # altitude - C = parmgi[alt_idx, :, :].T # extract data at fixed altitude - C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping - - print("Plotting each altitude cut...") - ax.plot_surface(X, Y, Z, facecolors=plt.cm.gnuplot2(C_normalized), rstride=1, cstride=1, edgecolor='none') - -print("Making combined altitude slice figure...") -# Add colorbar + labels -cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) -cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)', 'ne (m$^{-3}$)' -ax.set_xlabel('Latitude (deg)') -ax.set_ylabel('Longitude (deg)') -ax.set_zlabel('Altitude (m)') -ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want -ax.set_title('Electron Density at various altitudes') # change for whatever parm - -print("Saving figure...") -# Save + show plots -filename = 'ne_alt_cuts.png' # adjust filename for plot being saved -plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution - -plt.show() diff --git a/hz_processing_scripts/lat_cut_plotting.py b/hz_processing_scripts/lat_cut_plotting.py deleted file mode 100644 index 08fc688..0000000 --- a/hz_processing_scripts/lat_cut_plotting.py +++ /dev/null @@ -1,65 +0,0 @@ -import gemini3d.read as read -from gemini3d.grid.gridmodeldata import model2geogcoords -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.cm import ScalarMappable -from matplotlib.colors import Normalize - -# File setup -direc = '/Users/clevenger/Simulations/poker_asi_mar19/' # path to directory containing simulation output files -outdir = '/Users/clevenger/Simulations/poker_asi_mar19/postprocessing/' # path to save plot outputs - -# File reading stuff -cfg = read.config(direc) -xg = read.grid(direc) -parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 -dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index - -lalt = 128; llon = 128; llat = 128; # resolutions to plot model outputs (higher=better) - -print("Sampling in geographic coords...") -galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) - -print("Creating latitude slices...") -latitudes = [61, 63, 65, 67] # add/change latitudes as needed -lat_indices = [np.argmin(abs(glati - lat)) for lat in latitudes] - -# Lock longitude and altitude into place -ilon = np.argmin(abs(gloni - 212)) -ialt = np.argmin(abs(galti - 110)) - -# Creating 3D figure space to put 2D maps onto -fig = plt.figure(figsize=(10, 8)) -ax = fig.add_subplot(111, projection='3d') -ax.grid(False) - -# Normalize plot color/quantities so all of the 2D maps can share the same colorbar -norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) -mappable = ScalarMappable(norm=norm, cmap='viridis') - -print("Making latitude cuts...") -for lat_idx in lat_indices: - X, Z = np.meshgrid(gloni, galti) # longitude and altitude - Y = np.full_like(X, glati[lat_idx]) # latitude - C = parmgi[:, :, lat_idx].T # extract data at fixed latitude - C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping - - print("Plotting each latitude cut...") - ax.plot_surface(X, Y, Z, facecolors=plt.cm.viridis(C_normalized), rstride=1, cstride=1, edgecolor='none') - -print("Making combined altitude slice figure...") -# Add colorbar + labels -cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) -cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)' -ax.set_xlabel('Longitude (deg)') -ax.set_ylabel('Latitude (deg)') -ax.set_zlabel('Altitude (m)') -ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want -ax.set_title('Electron Density at various latitudes') # change for whatever parm - -print("Saving figure...") -# Save + show plots -filename = 'ne_lat_cuts.png' # adjust filename for plot being saved -plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution - -plt.show() diff --git a/hz_processing_scripts/lon_cut_plotting.py b/hz_processing_scripts/lon_cut_plotting.py deleted file mode 100644 index 9b67496..0000000 --- a/hz_processing_scripts/lon_cut_plotting.py +++ /dev/null @@ -1,65 +0,0 @@ -import gemini3d.read as read -from gemini3d.grid.gridmodeldata import model2geogcoords -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.cm import ScalarMappable -from matplotlib.colors import Normalize - -# File setup -direc = '/Users/clevenger/Simulations/poker_asi_mar19/' # path to directory containing simulation output files -outdir = '/Users/clevenger/Simulations/poker_asi_mar19/postprocessing/' # path to save plot outputs - -# File reading stuff -cfg = read.config(direc) -xg = read.grid(direc) -parm = "ne" # choose ne, Ti, Te, v1, v2, v2, J1, J2, or J3 -dat = read.frame(direc, cfg["time"][2], var=parm) # change [2] for desire index - -lalt = 128; llon = 128; llat = 128; # resolutions to plot model outputs (higher=better) - -print("Sampling in geographic coords...") -galti, gloni, glati, parmgi = model2geogcoords(xg, dat[parm], lalt, llon, llat, wraplon=True) - -print("Creating longitude slices...") -longitudes = [200, 205, 210, 215, 220] # add/change latitudes as needed -lon_indices = [np.argmin(abs(gloni - lon)) for lon in longitudes] - -# Lock latitude and altitude into place -ilat = np.argmin(abs(glati - 65)) -ialt = np.argmin(abs(galti - 110)) - -# Creating 3D figure space to put 2D maps onto -fig = plt.figure(figsize=(10, 8)) -ax = fig.add_subplot(111, projection='3d') -ax.grid(False) - -# Normalize plot color/quantities so all of the 2D maps can share the same colorbar -norm = Normalize(vmin=np.nanmin(parmgi), vmax=np.nanmax(parmgi)) -mappable = ScalarMappable(norm=norm, cmap='viridis') - -print("Making longitude cuts...") -for lon_idx in lon_indices: - X, Z = np.meshgrid(glati, galti) # latitude and altitude - Y = np.full_like(X, gloni[lon_idx]) # longitude - C = parmgi[:, lon_idx, :].T # extract data at fixed longitude - C_normalized = (C - np.nanmin(C)) / (np.nanmax(C) - np.nanmin(C)) # normalize for color mapping - - print("Plotting each longitude cut...") - ax.plot_surface(X, Y, Z, facecolors=plt.cm.viridis(C_normalized), rstride=1, cstride=1, edgecolor='none') - -print("Making combined longitude slice figure...") -# Add colorbar + labels -cbar = fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=10) -cbar.set_label('ne (m$^{-3}$)') # or 'Ti (K)', 'Te (K)', 'v$_{1} (m/s)$', 'J_{1} (A/m$^{2}$)' -ax.set_xlabel('Latitude (deg)') -ax.set_ylabel('Longitude (deg)') -ax.set_zlabel('Altitude (m)') -ax.view_init(20, 45) # adjust to move 3d space around to view cuts however you want -ax.set_title('Electron Density at various latitudes') # change for whatever parm - -print("Saving figure...") -# Save + show plots -filename = 'ne_lon_cuts.png' # adjust filename for plot being saved -plt.savefig(f'{outdir}/{filename}', dpi=300, bbox_inches='tight') # adjust dpi for plot resolution - -plt.show() \ No newline at end of file