import xarray as xr
import arcpy
from pathlib import Path
from datetime import datetime, timeWorkflow for STAQS GCAS netCDF files in ArcGIS Pro
Summary
This notebook shows how ArcGIS users can access Synergistic TEMPO Air Quality Science (STAQS) data using ArcGIS PRO.
When accessing data files (in Network Common Data Form (netCDF) format) from the STAQS field campaign in ArcGIS Pro, users have reported issues with improper display of data from the Langley Research Center’s (LaRC’s) GIII GeoCAPE Airborne Simulator (GCAS) instrument. The goal of this notebook is to help with properly displaying the data in ArcGIS Pro. The workflow shown in this notebook creates new netCDF files with the desired data while removing null values and displaying the new netCDF files in ArcGIS Pro. More information about each step is available in the comments for each cell.
Prerequisites
- ArcGIS Pro
arcpyxarray
Setup
Import required packages. If you don’t have them, run “pip install {package_name}” in your terminal.
User needs to specify the following paths
- Assign directories and variables, you will need to change your projects Geodatabase path, keep it in the same format as mine
- For the sake of this workflow, add the NC_Files and New_NC_Files directories to your desired projects folder
nc_dir = "NC_Files/"
output_dir = "New_NC_Files/"
geo_database = r"C:\Users\gmojica\Documents\ArcGIS\Projects\RasterTestMosaic\RasterTestMosaic.gdb"input_path = Path(nc_dir)
output_path = Path(output_dir)
output_list = list(output_path.iterdir())aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.listMaps()[0]Iterate through directories to modify nc files and add them to the New_NC_Files dir, assign variables, change what you need to.
lat = "lat"
lon = "lon"
var = "no2_vertical_column_below_aircraft"
t = "time"
xt = "xtrack"
for file in list(input_path.iterdir()):
dataset = xr.open_dataset(file)
dataset = dataset[[var, lat, lon]].dropna(dim=t)
file_name = str(file).replace(str(input_path), "")
new_dataset = dataset.to_netcdf(path=str(output_path) + str(file_name))
aprx.save()
print(f"New NetCDF created at {datetime.now()}")New NetCDF created at 2025-05-12 11:53:19.441491Run Geoprocessing tool to create NetCDF Feature Layer
- The NetCDF Feature Layer is a temp layer and cannot be accessed, need to use Copy Feature tool to create a permanent feature layer
- If this cell doesn’t produce a feature layer, you may need to either refresh the folder with the new nc file and wait a few minutes, or save the project, close, and reopen ArcGIS Pro and run the first few cells again.
- This cell has experienced issues before with freezing up, if it takes longer than a few minutes to run (depending on the amount of files), you may want to perform a quick restart of your system, that should free up any bottlenecks with this cell.
print("Creating NetCDF Feature Layers, this may take some time.")
for file in output_list:
arcpy.md.MakeNetCDFFeatureLayer(
in_netCDF_file=str(file),
variable=var,
x_variable=lon,
y_variable=lat,
out_feature_layer=str(file) + "_NC_Flyr",
row_dimension=f"{t};{xt}",
z_variable="",
m_variable="",
dimension_values=None,
value_selection_method="BY_VALUE",
)
aprx.save()
print(f"NetCDF Feature Layer(s) created at {datetime.now()}")Creating NetCDF Feature Layers, this may take some time.
NetCDF Feature Layer(s) created at 2025-05-12 11:53:56.192013Rename temp feature layer to allow for Copy Features tool to run successfully
Note: Once you run this cell once, do not run it again.
lyr_list = []
for layer in m.listLayers("New_NC*"):
layer.name = str(layer).replace("New_NC_Files\\", "").replace(".nc", "").replace("-", "")
lyr_list.append(layer.name)
aprx.save()
print("Layer renamed, you can create permanent feature class.")Layer renamed, you can create permanent feature class.Copy Feature to create permanent feature classes.
print("Creating permanent features, this make take some time.")
for layer in m.listLayers("staqs*"):
arcpy.management.CopyFeatures(
in_features=layer,
out_feature_class=str(geo_database) + "\\" + str(layer),
config_keyword="",
spatial_grid_1=None,
spatial_grid_2=None,
spatial_grid_3=None,
)
arcpy.Delete_management(layer)
print(f"Done at {datetime.now()}")Creating permanent features, this make take some time.
Done at 2025-05-12 11:54:35.129398Use search cursor to create datetime.date variable and then run select analysis.
The time variables need to be input in this format: (HH, MM, SS). If it’s a single digit, do not precede with 0, only use the single digit (1, not 01).
feature_classes = arcpy.ListFeatureClasses()
fc_field = t
start_time = time(15, 30, 0)
end_time = time(18, 0, 0)
print(
"Creating feature class(es) based on time bounds, deleting input features. This may take some time."
)
for fc in feature_classes:
with arcpy.da.SearchCursor(fc, fc_field) as cursor:
for row in cursor:
date = row[0].date()
break
upper_bound = datetime.combine(date, start_time)
lower_bound = datetime.combine(date, end_time)
input_feature = fc
output_feature = str(geo_database) + "\\" + str(fc) + "_copy"
where_clause = f"{t} >= timestamp '{upper_bound}' AND {t} <= timestamp '{lower_bound}'"
arcpy.analysis.Select(input_feature, output_feature, where_clause)
arcpy.Delete_management(fc)
aprx.save()
print(f"Done at {datetime.now()}")Creating feature class(es) based on time bounds, deleting input features. This may take some time.
Done at 2025-05-12 11:54:50.864299Interpolate by natural neighbor with z value of your variable, feel free to change cell size.
cell = 0.00707232999999999
print("Interpolating by Natural Neighbor.")
for layer in m.listLayers("staqs*"):
arcpy.ddd.NaturalNeighbor(
in_point_features=str(layer),
z_field=str(var),
out_raster=str(geo_database) + "\\" + str(layer) + "_intNN",
cell_size=cell,
)
aprx.save()
print(f"Done at {datetime.now()}")Interpolating by Natural Neighbor.
Done at 2025-05-12 11:54:58.926297Run the following cell if you’d like to delete the NetCDF Feature Layer and keep the interpolated raster layer
If you’d like to keep the Feature Layers on the map, don’t run this cell. If you delete the feature layer and want it back, it is stored in your projects GDB.
for m in arcpy.ListFeatureClasses():
arcpy.Delete_management(m)
aprx.save()
print("Feature layer(s) have been removed from the map.")Feature layer(s) have been removed from the map.