Where is normalisation by cosine latitude

[cjmerchant]

The code is complex, so maybe I miss it but for calculating the latitudinal density histogram I see only a line like

n_sst_q4_lat = np.sum(ds['n_sst_q4_lat'],axis=0)[0:3600,] / np.array((lat_fraction * ocean_area) / years)

lat_fraction appears to be the fraction of the area of that latitude bin that is sea

however, the equation is only correct is lat_fraction is the fraction of the total ocean area that is in the latitude range, which is not the same thing

[patternizer]

lat_fraction is returned from the method: calc_lat_fraction(). It calculates the fraction of non-land pixels by applying the water, sea-ice, land mask from the OSTIA L4 product from the ESA SST CCI project, produced using OSTIA reanalysis sytem v3.0:

The fraction of non-land pixels (i.e. water + sea ice) in each latitudinal slice is then calculated as follows (with x=ds.lon):

f = 1 - (np.sum(land[0,:,:],axis=1) / len(x)*1.)

Since the resolution of lat_vec is the same as the NAVOCEANO_landmask_v1.0 EUMETSAT_OSI-SAF_icemask ARCLake_lakemask, no cosine normalisation is performed and no interpolation is performed.

[cjmerchant]

So, this is as I thought. For each latitude zone "ocean_area", which is a single constant, needs to be replaced with "total_area_of_surface_in_this_latitude_zone" which is approximately

length width = (2piradius_of_Earthcosine(latitude_at_centre_of_zone)(piradius/latitude_width_of_zone)

As a check, we should be able to find the weighted mean of the latitude zone plot and check it is the same as the mean density found by looking at the total sume over the whole ocean area:

w_mean = (SUM of (zonal_density (cosine(latitude_at_centre_of_zone)lat_fraction))) / (SUM of (cosine(latitude_at_centre_of_zone)*lat_fraction)))

[cjmerchant]

I believe np.array((lat_fraction ocean_area) / years) is wrong in another way too. We want the density = number per area per year = km-2 year-1 But the above gives units of km-2 year So it should be np.array((lat_fraction ocean_area) * years)

[patternizer]

I believe np.array((lat_fraction ocean_area) / years) is wrong in another way too. We want the density = number per area per year = km-2 year-1 But the above gives units of km-2 year So it should be np.array((lat_fraction ocean_area) * years)

I'm running it through the de-bugger now and doing some checks (thanks). You are correct, thank you.

[patternizer]

So, this is as I thought. For each latitude zone "ocean_area", which is a single constant, needs to be replaced with "total_area_of_surface_in_this_latitude_zone" which is approximately

length width = (2_pi_radius_of_Earth_cosine(latitude_at_centre_ofzone)(piradius/latitude_width_of_zone)

As a check, we should be able to find the weighted mean of the latitude zone plot and check it is the same as the mean density found by looking at the total sume over the whole ocean area:

w_mean = (SUM of (zonal_density (cosine(latitude_at_centre_of_zone)lat_fraction))) / (SUM of (cosine(latitude_at_centre_of_zone)*lat_fraction)))

The equation for the area of the Earth between a line of latitude and the north pole (the area of a spherical cap) is:

A = 2piR*h

where R is the radius of the earth and h is the perpendicular distance from the plane containing the line of latitude to the pole. We calculate h using trigonometry:

h = R*(1-sin(lat))

The area north of a line of latitude is therefore:

A = 2piR^2(1-sin(lat))

The area between two lines of latitude is the difference between the area north of one latitude and the area north of the other latitude, i.e.:

A = | 2piR^2(1-sin(lat2)) - 2piR^2(1-sin(lat1)) = 2piR^2 |sin(lat1) - sin(lat2)

Our lat_vec is: [-89.75, 89.57, step=0.05], therefore, to perform the above calculation centered on these values, we need to obtain the sin(lat) at lat +/- (0.05/2). This has been added into the code now.

This is how the ocean area scales within the cosine latitude envelope:

The resulting histogram of SST(lat) is: