how to show the number of grid cells and calculate area

[pybl]

Hi,

I am trying to compute the surface area of a selected region, or at least knowing the number of grid cells so that I can compute surface area by myself.

for example, if using xarray:

dataset.TS.where(Temperature >273).plot() will show me the area where surface temperature is greater than freezing point.

But, how can I get to know the number of grid cells?

Thanks in advance for your help!

 

[mlevy]

If I understand your question correctly, you want to do something like

>>> np.sum(np.logical_not(np.isnan(ds.TS.where(ds.TS > 273))))
<xarray.DataArray 'TS' ()>
array(34162)

In this example, there are 192*288=55296 grid cells and 34162 of them have a temperature above 273 K. I don't see how this helps you calculate the surface area of that region, though, because the cells are not all the same size. I'm going to move this to the CAM forum and see if someone there can help you find the area in the history file... at that point, the area you want would be

>>> area_da.where(ds.TS > 273).sum()

 

[pybl]

Thank you so much for the reply!

The resolution of every cell is the same, the area depends on their latitude, if I understand it correctly, due to the curvature of Earth

I am thinking of applying gaussian weights to the grid cells at different latitudes to calculate the area of the grid cells.

If I can do something like this (the format may be wrong):

step 1: filter out TS = ds. TS.where(ds.TS > 273)
step 2: define resolution in degrees: 'res' = horizontal resolution * vertical resolution (should be a constant about ~ 1 degree squared)
step 3: assign the resolution value to each grid in TS by replacing the surface temperature value with 'res', let's call it 'res_map'

Now, I would have a data_array with a coordinate map that relates to TS > 273 K and each grid cell has the value of 'res'

step 4: define the function 'guassian_weight' and apply it to 'res_map', to get the scaled resolution for each cell, let's call it 'scaled_res_map'
step 5: sum up scaled_res_map.values, and multiply by the conversion between the area in degree squared and the area in kilometer squared

Please point out:
- if there is a better solution than this because I think CAM has a parameter called 'area' but I did not find it in my history file.
- or anywhere you think the above-proposed method is wrong.

Thank you for helping!

 

[cacraig]

I spoke with one of our scientists and here is their response:

I usually grab the areas from the domain files. They are also default i/o for se-dycore on the native grid, but not when you have interpolate_output=T. Although for lat-lon grid you can compute area using default cam i/o "gw," the guassian weights

An example of computing area on lat-lon grids is in example 1 here: https://www.ncl.ucar.edu/Document/Functions/Built-in/wgt_areaave.shtml

 

[pybl]

Thanks a lot! :)