Face area calculations for sphere radius > 1

[st-ncar]

I work with MPAS data on a sphere of radius $r_E > 1$ and look for ways to compute integrals of quantities defined at the cell centers of the primal mesh. This requires computing the face areas (either explicitly or implicitly as part of UxDataArray.integrate()) and I am not sure I fully understand the behavior of uxarray. There seem to be at least two ways to obtain the face areas:

• Grid.compute_face_areas(): computes the face area from the node coordinates (spherical or Cartesian).
• Grid.face_areas: returns face_areas if contained in Grid, otherwise calls the function above

Intuitively, I would expect that both methods (with either coordinate system) give the same answer. When I tested this, I got a different answer in each case:

1. Grid.compute_face_areas(latlon=True) seems to give the face areas on the unit sphere ($r = 1$)
2. Grid.compute_face_areas(latlon=False) gives very small face areas (1e-16 or smaller). I wasn't sure how to interpret this.
3. Grid.face_areas seems to give the face area on the sphere with radius $r = r_E$ (perhaps these areas are read from the grid file directly?)

I can rescale the result of method 1 to recover method 3. It also seems that UxDataArray.integrate() uses method 1 internally and I have to rescale the result to get the expected integral value on a sphere with radius $r_E$.

I was wondering if someone can help me interpret the result of the 3 methods above (what I wrote are just my guesses). Also, do you have a recommendation for which method I should use? And in case the above interpretation is correct, would it be possible to make Grid.compute_face_areas() aware of Grid.attrs['sphere_radius'] so that there is no need for rescaling between method 1 and 3?

[philipc2]

@rajeeja

[rajeeja]

Yes, will update soon

[st-ncar]

Great thanks, looking forward to the next release 😄 uxarray is a great tool, thanks for all your work!

[philipc2]

We should consider adding a Grid.sphere_radius property that can be set by the user and/or parsed from the input grids.

This would allow us to tweak the default behavior of Grid.bounds to use that internal property.

[rajeeja]

I was wondering if someone can help me interpret the result of the 3 methods above (what I wrote are just my guesses). Also, do you have a recommendation for which method I should use? And in case the above interpretation is correct, would it be possible to make Grid.compute_face_areas() aware of Grid.attrs['sphere_radius'] so that there is no need for rescaling between method 1 and 3?

Do you have a small file to reproduce this?

[rajeeja]

I work with MPAS data on a sphere of radius r E > 1 and look for ways to compute integrals of quantities defined at the cell centers of the primal mesh. This requires computing the face areas (either explicitly or implicitly as part of UxDataArray.integrate()) and I am not sure I fully understand the behavior of uxarray. There seem to be at least two ways to obtain the face areas:

• Grid.compute_face_areas(): computes the face area from the node coordinates (spherical or Cartesian).
• Grid.face_areas: returns face_areas if contained in Grid, otherwise calls the function above

Intuitively, I would expect that both methods (with either coordinate system) give the same answer. When I tested this, I got a different answer in each case:

1. Grid.compute_face_areas(latlon=True) seems to give the face areas on the unit sphere (
   r

   1
   )

2. Grid.compute_face_areas(latlon=False) gives very small face areas (1e-16 or smaller). I wasn't sure how to interpret this.

3. Grid.face_areas seems to give the face area on the sphere with radius
   r

   r
   E
   

   (perhaps these areas are read from the grid file directly?)

I can rescale the result of method 1 to recover method 3. It also seems that UxDataArray.integrate() uses method 1 internally and I have to rescale the result to get the expected integral value on a sphere with radius r E .

I was wondering if someone can help me interpret the result of the 3 methods above (what I wrote are just my guesses). Also, do you have a recommendation for which method I should use? And in case the above interpretation is correct, would it be possible to make Grid.compute_face_areas() aware of Grid.attrs['sphere_radius'] so that there is no need for rescaling between method 1 and 3?

my gut feeling is one assumes that the sphere is unit radius and the other doesn't; we are in the process of removing the latlon variable, the implementation will be in cart coords and conversion will happen as a part of it internally.

[st-ncar]

Thanks for looking into this! Do you have access to the Glad file service? If so, you can find an example data set in /glade/derecho/scratch/stoedtli/test_psfc/60km_uniform/from_forecast/. The relevant files in this directory are:

• test_integral.ipynb: short notebook that illustrates the different results for the 3 methods
• invariant.163842.nc: file that contains the example MPAS mesh
• p_sfc.*.nc: an example data set defined over the MPAS mesh (not strictly needed for the surface area calculations. I became aware of the differences in surface area calculations when trying to integrate quantities in this data set)

[philipc2]

Thanks for sharing this! I do have access.

[philipc2]

I can see the notebook and data set, however I can't access the invariant.163842.nc grid file.

Could you chmod o+r the file? This should give me read permission.

[st-ncar]

Sorry about that! The original invariant.163842.nc file was a symbolic link, I have now replaced it with the actual file. I have also updated the file permissions. Please let me know if you still can't access it.

[philipc2]

It's good now! Thank you!