ATS observations at multiple points

Sundar Niroula

unread,

Apr 26, 2024, 12:39:07 PM4/26/24

to Amanzi-ATS Users

Hi ATS community,

I am trying to extract some ATS simulated variables as observations say, transpiration and permeability, at hundreds of points (after identifying the horizontal coordinates of the surface) as a csv file. I know this can be extracted from ats_vis_data, but for now I am looking for other options.

Second approach is possibly by defining a region as “point” for a coordinate in the XML and repeating for each of them. But doing this for hundreds of points will overwhelm the xml file as I need the observations for more than two variables eventually.

I wonder if there is a third approach, which is possibly via adding these hundreds of points in the mesh as labeled set in WW. I followed a similar discussion regarding this. I tried to adapt this approach but could not succeed. Here are my steps in WW:

1. First, the lat and lon of points of my interest are stored in df_NHD_points. There are 528 observation points. m2 is my mesh.

2. Appended the observation_points as cell in the mesh.

all_observation_point_list = []
for index, row in df_NHD_points.iterrows():
# Create a shapely Point object from Longitude and Latitude coordinates
point = Point(row['Longitude_RC'], row['Latitude_RC'])
# Append the Point object to the list
all_observation_point_list.append(point)

# generate the corresponding observation cells
my_obs_cell=[]
# Searching for the index
for ic, cell in enumerate(m2.conn):
points = np.array([m2.coords[i] for i in cell])
polygon = shapely.geometry.Polygon(points)
# Check if any point in all_observation_point_list is contained within the polygon
for observation_point in all_observation_point_list:
if polygon.contains(observation_point):
my_obs_cell.append(ic)

3. Added the labeled set and extruded the mesh

# add the observation cells to 2D mesh
new_id = m2.next_available_labeled_setid()
m2.labeled_sets.append(LabeledSet("obs_cells", new_id, "CELL", my_obs_cell,
to_extrude=True))

Doing this, I can see the labeled set and ent_ids allocated for “obs_cell” for both 2D and 3D. “obs_cell” also appears in XML as region: labeled set.

I first tried to get “transpiration” as follows:

The simulation went through, generated “obs_cell_transpiration.csv”. But all values for transpiration were “nan”. Also, although I got all nan values at daily time step, I got only 1 vector/column for transpiration. But I was expecting 528 columns, as I have 528 points defined in obs_cells. In other words, the output matrix was 1 X Time while I am looking for 528 X Time.

My Questions:

1) I might be missing many things here. Can you please help me identify? More importantly, is this (the third) approach reasonable?

2) I assume I need to define a new .csv file to store the output for every new variable of my interest. Is this assumption valid?

2) Can this (third) approach be helpful if I need the outputs at these observation points for their entire column, e.g., porosity?

Any suggestions will be helpful.

Thank you,

Sundar Niroula

Amanzi-ATS Users

unread,

Apr 26, 2024, 5:53:09 PM4/26/24

to Amanzi-ATS Users

I think your first try (using region: point) was the right way to go -- you need to define a region for every point, and an observation for every variable at every region. You would have to do this no matter how you defined the region (as a point or as a labeled set). Putting all the points in a single labeled set just means you'll be averaging (or otherwise integrating) the values across all points and saving one observation.

These can all be put in the same csv if you want, but I'd probably keep each variable as its own file -- 528 columns is enough, no need for NVARIABLES * 528 columns.

I would suggest scripting the generation of the xml for both regions and observations. Using amanzi_xml (from $AMANZI_SRC_DIR/tools/amanzi_xml) you can add an observation file entry, then loop over each region and variable and add the column.

The attached script is (untested) pseudocode that should give you a hint of how to proceed with this.

Ethan

observations.py

Message has been deleted

Sundar Niroula

unread,

Apr 29, 2024, 7:06:23 PM4/29/24

to Amanzi-ATS Users

Hi Ethan,

Thank you for your response. I adapted the script you shared and it worked fine.

Thanks,

Sundar Niroula

unread,

May 6, 2024, 7:16:44 PM5/6/24

to Amanzi-ATS Users

Hi all,

Following up on this thread,

While I was able to accurately extract surface observations by defining {x,y} as region point (with Ethan's help, Thanks!), I am getting “nan” when I extract subsurface variables like “water content” and “saturation-liquid’ across depths. The nan values are only seen in the observation (.csv) files. The model runs successfully, and I can see non-nan values for “water content” and “saturation-liquid’ in paraview.

For subsurface observations, I use {x,y,z} as region point. I define z as the centroid elevation of the cell that contains {x,y}, given by z = m2.centroid[xy_containing_cell_index,2]. For layer 0, I simply assign z as the elevation (to see what happens!). For the first layer (layer1), I set elevation (or z1) to be the difference between z (surface) and half the depth of the first layer. (Just to note, my first layer is 0.05m thick). And so on for other layers. The following is an example of my approach for layer 0 and layer 1 at one point location.

Region Definition:

Observation Definition:

</ParameterList>

</ParameterList>

<ParameterList

Doing so, generates a .csv file with nan at all time period, all points, and all columns (mesh) layers. The result is same for “saturation-liquid” too. I also tried changing functional = “average” but no change.

I might be missing something here or making inaccurate assumptions. Has anybody faced and resolved such issue? I know there was a similar discussion here before but there were no further updates.

I would appreciate any thoughts or suggestions.

Thank you,
Sundar Niroula

Coon, Ethan

unread,

May 7, 2024, 1:01:32 PM5/7/24

to Sundar Niroula, Amanzi-ATS Users

So nan means that the region was empty (as it divides by total cell volume for “integral” or “average”), which it sounds like you concluded too.

Your process sounds right for getting z coordinates that should be “in the prism” defined in the 3D mesh. We’ve had some trouble in the past of getting the cell’s “does a cell contain a given point” algorithm correct, so it’s possible you’re doing it right and the code is still not resolving the point to the right cell. But I thought that was fixed. Can you plot the mesh in Paraview or VisIt and confirm that the point you are using for the region is correct (e.g. “inside the mesh and should resolve to _some_ cell”). If you can do that, send me directly a mesh and input file, and I can get someone to take a look to see why the point region isn’t being found in any cell.

Ethan

--
You received this message because you are subscribed to the Google Groups "Amanzi-ATS Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ats-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ats-users/f0b8a8c1-2694-4bc7-8dc1-b59940c49346n%40googlegroups.com.

Wyatt Tatge

unread,

Mar 9, 2026, 5:02:32 PM (4 days ago) Mar 9

to Amanzi-ATS Users

Hi all,

I am attempting to output subsurface pressure at multiple points. At my field site, I have 12 wells that are multilevel, that are screened at different intervals, each with unique water level data. I got the elevations for each well from watershed workflow by using z = m2.centroids[obs_wells,2] (similar to above), and I double checked these elevations by writing an output from ATS and manually in Paraview. I input each well screen as a point in ATS as seen below:

<ParameterList name="w6p5" type="ParameterList">
<Parameter name="region type" type="string" value="point"/>
<Parameter name="coordinate" type="Array(double)" value="{253874.0, 4927916.0,320.03}"/>
</ParameterList>
<ParameterList name="w6p6" type="ParameterList">
<Parameter name="region type" type="string" value="point"/>
<Parameter name="coordinate" type="Array(double)" value="{253874.0, 4927916.0,318.51}"/>
</ParameterList>
<ParameterList name="w6p7" type="ParameterList">
<Parameter name="region type" type="string" value="point"/>
<Parameter name="coordinate" type="Array(double)" value="{253874.0, 4927916.0,317.69}"/>
</ParameterList>

However, for each of the wells in the above example, my ATS output has the same pressure. I also output depth from ATS, and each shows a depth of 0.0,5 putting them all in my first layer cell. A similar thing is happening with a few of my other wells. The well in the example had a surface elevation of 321.86. I have also used Paraview to check the visualization at each of these points, and they correspond to different cells with different pressures. It seems to me that ATS is not correctly identifying the location of these points. I'm not sure if anyone has experienced this or has an idea of how to fix it. I am using ATS-1.6. Here is how I am outputting each of these in the observation as well:

<ParameterList name="w6p5 pressure [Pa]" type="ParameterList">
<Parameter name="variable" type="string" value="pressure"/>
<Parameter name="region" type="string" value="w6p5"/>

Thanks
Wyatt

Bo Gao

unread,

Mar 9, 2026, 10:51:19 PM (3 days ago) Mar 9

to Wyatt Tatge, Amanzi-ATS Users

Hi Wyatt,

Can you please attach your mesh, input file, and output observation file? I tested v1.6 without problems for point observations.

Best,

Bo Gao

To view this discussion visit https://groups.google.com/d/msgid/ats-users/38ce16ef-4385-4dfb-ac6c-80544f70f499n%40googlegroups.com.

Wyatt Tatge

unread,

Mar 10, 2026, 9:24:09 AM (3 days ago) Mar 10

to Bo Gao, Amanzi-ATS Users

Hi Bo,

I have attached the mesh, input file, and outputs in the zip folder.

Thanks

Wyatt

sim_results.zip

Bo Gao

unread,

Mar 11, 2026, 6:46:52 PM (2 days ago) Mar 11

to Wyatt Tatge, Amanzi-ATS Users

Hi Wyatt,

I guess there should be some issues to locate the cell by a point coordinate. For the w6p5, w6p6, w6p7 three points in your example, all return the top cell id. I do not know how point locating works in the code. For your case, I would recommend using a very small box instead of a point to define a region for a well's point, and use average for functional in the observation.

Best,

Bo Gao

Coon, Ethan

unread,

Mar 11, 2026, 7:06:39 PM (2 days ago) Mar 11

to Bo Gao, Wyatt Tatge, Amanzi-ATS Users

Hmm, I thought we fixed these issues. I’ll see what I can find.

Ethan

From: ats-...@googlegroups.com <ats-...@googlegroups.com> on behalf of Bo Gao <boga...@gmail.com>
Date: Wednesday, March 11, 2026 at 4:47 PM
To: Wyatt Tatge <wyatt...@yale.edu>
Cc: Amanzi-ATS Users <ats-...@googlegroups.com>
Subject: Re: [EXTERNAL] Re: ATS observations at multiple points

This Message Is From an External Sender

This email was sent from a non-ORNL address. If suspicious, use the Report Phish button in Outlook.

To view this discussion visit https://groups.google.com/d/msgid/ats-users/CAN8t3fK%3D3EM-_%3D9H%3D35WggzYzJE5fGeThR20KgoTzi%3DXqGEONA%40mail.gmail.com.

Reply all

Reply to author

Forward