upper and lower area calculated from mask do not equal

[miniufo]

This is raised by @AminFazlKazemi.

well If the equivalent latitude is chosen perfect, you can fill in values painted as light blue in white space below Equivalent Latitude so The resultant Adiabaticcaly sorted filed is produced. I chose a mask output from your code corresponding e.g. Latitude=45 degrees North. Pixels are reclassied there with -1 ,+1 and zero. I then multiplied each pixel by its corresponding fractional area in squared kilometers. If The choosing of Equivalent Latitude is perfect, then integrating area over pixels shown with +1 is exactly equals area over pixels shown with -1. Of course the resolution of dataset plays a role here.

I did not check this carefully. Maybe it is related to the discretized grid points or weights calculated in xgcm. But it needs to be double check.

I made a rough estimate as follows. Grid resolution is 0.7016754 degree. So the area of a grid point at 45°N is (0.7016754*111.11)**2*cos(45°), which is about 4297.32km^2. Similarly, area is 3039.13km^2 at 60N and 1029.35km^2 at 80.25N. It is reasonble to take this as the grid-scale error. I guess your unit of km^2 is wrong (should be m^2 for 10^12).

Also, how do you calculate the area? can you post your code snippet here? Maybe the algorithms are also somewhat different.

[AminFazlKazemi]

Hi

You are correct the units are squared meters.

Area=(numpy.radians(resolution)*earth_radius)**2 *(numpy.cos(numpy.radians(data.latitude)))+0*data.longitude

well , yes I used earth-radius as not complete sphere but the effect of this is partial. I checked code of hn2016_falwa packages. It works just fine even near poles of course between 85S to 85 N (it is confined in fortran codes as I had done in my own code) but of course not exactly.

[miniufo]

Here is what I thought. The (normalized) area bounded by neighbouring lats and lons is exactly defined as: |lon1 - lon2| * |sin(lat1) - sin(lat2)| / (4*PI) See here for a schematic illustration. We can call this method 1.

However, xcontour is built on xgcm which uses discretized way of estimating area of each grid point. Which should be close to the formula you used (dx * dy) * cos(lat). We can call this method 2.

So you could try the following things:

1. Interpolate the tracer onto a much higher resolution using xarray's interp (given finer coordinates). Then compare the areas of W+ and W- from masks (after interpolation) and see if the differences would reduce when resolution is higher;
2. Use the internal method of xgcm to estimate areas of W+ and W- as areaWP = grid.integrate(maskWP, dims=['X', 'Y']) and areaWN = grid.integrate(maskWN, dims=['X', 'Y']). Then see if the differences are reduced.
3. You can also check the way of calculating area in Huang's Fortran program and see if it is method 1.

I do not get your data and hope you can help me find the reasons and improve the accuracy of this package. Thank you very much.

[miniufo]

I just tried with my case and Huang's barotropic case, confirming that if this line:

m = eqDim>eqDim.values[j] if coord_incre else eqDim<eqDim.values[j]

is changed to:

m = eqDim>=eqDim.values[j] if coord_incre else eqDim<=eqDim.values[j]

then the resulted areas of W+ and W- are much closer. Please make this small change and see if your case improves.

[AminFazlKazemi]

Hi. for Latitude=75 degrees North negative_intrusion_area= 9.08332754e+12 squared meters positive_intrusion_area = 4.88416361e+11 squared meters

[AminFazlKazemi]

This is the code I used. Itried but didn't manage to add indentations correctly. The data can be downloaded from era interim dataset with 0.75 degrees resolution for April 2011 on isentropic vertical levels. I targeted 320 kelvin but you can change it

[AminFazlKazemi]

LWA.txt This is the code

[AminFazlKazemi]

In the finite amplitude limit, we see enclosed or cut-off contours forming. I suggest test definition of area within contours visually in such a case. i also suggest to read Ghinassi2018 paper page 4103 part 3

In fact you need to create a polygon based on contour intervals and calculate area of each polygon. In GIS terminology, it is referred to as zonallization , e.g. "zonal statistics" in QGIS/Esri Arc GIS softwares. The polygon can be used with rasterization to define polygon number for each grid point. I had found a code that did this job. If I found it I will send. maybe GIS packages in python can help in this regard too. This will definitely solve the issue of grid resolution, some thing that grid based area summing-up cannot solve . The more exactly the the contours are defined, the more exact results are reached. Contour integration is one of the most advanced and powerful methods in meteorology and Earth sciences. I hope you can solve the issue so we cqan use your functiuons in real applications.

[miniufo]

Can you share you netcdf data so that I can reproduce your results?

[miniufo]

I've downloaded it and I may need some time to look into it. You can remove the link now.