Re: area_global_rectilinear_grid

Actually, there is an undocumented function in shea_util.ncl

;-------------------------------
undef("area_polar_region")
function area_polar_region (lat[1]:numeric, opt[1]:logical)
; Area from pole to latitude.
local rearth, pi, rad, rr, twopi, parea
begin
   rearth = 6371.220d0 ; default; kilometers
   if (opt .and. isatt(opt, "rearth")) then
       rearth = opt@rearth
   end if

   pi = 4d0*atan(1.0d0)
   rad = pi/180d0
   twopi = 2.0d0*pi
   parea = twopi*rearth*rearth*(1.0d0-sin(rad*lat))
   parea@long_name = "Polar Area"
   parea@units = "km^2"

   return(parea)
end

On 10/8/13 5:49 AM, Dennis Shea wrote:
> This is *NOT* correct. Not sure what I was thinking.
>
> Sorry
> D
>
> On 10/7/13 11:47 AM, Dennis Shea wrote:
>> As noted the function is an approximation.
>>
>> The surface area of a sphere is
>>
>> A = 4*pi*R^2 ; area of sphere; R is radius of earth
>>
>> Area of a hemisphere:
>>
>> diff = 90-0.0
>> Ahem = A*(diff/180.0) = A*0.5
>>
>> Area of the 'cap'
>>
>> diff = 90-88.92773
>> A889 = A*(diff/180.0)
>>
>> On 10/2/13 7:37 PM, Soumik Basu wrote:
>>> Hi,
>>>
>>> I am trying to calculate the sea ice extent from sea ice fraction data
>>> for CAM3.1 in T85 grid.
>>> I used the area_global_rectilinear_grid function to calculate the area
>>> of each of the grid box. But I am trying to estimate the area over the
>>> pole as the latitudes ends at 88.92773. For doing that I tried
>>> calculating in the following way. But what confused me is that why the
>>> sum of the grid areas calculated using the function is greater than the
>>> surface area of the earth. Is there any better way to estimate the area
>>> of the small circle over the pole?
>>>
>>>
>>> ;****************************************************
>>> ; Calculate area at each grid point
>>> ;****************************************************
>>> ; T85
>>>
>>> nlat = 128
>>> lat = latGau(nlat, "lat", "latitude", "degrees_north")
>>> mlon = 256
>>> lon = lonGlobeF(mlon, "lon", "longitude", "degrees_east")
>>> print(lat)
>>> area = area_global_rectilinear_grid (lat, lon, False)
>>>
>>> areaa = dble2flt(area)
>>>
>>> printVarSummary(areaa)
>>> printMinMax(areaa, True)
>>>
>>> R = 6371.220 ; in km
>>>
>>> pi = 4d0*atan(1d0)
>>>
>>> S_Area = 4d0*pi*R^2 ;Total surface area of the earth
>>>
>>> print(S_Area)
>>>
>>> Area_total = sum(areaa) ; Total area from the grid boxes
>>>
>>> print(Area_total)
>>>
>>> Pole_area = (S_Area - Area_total)/(10^6)*2 ; Area over the pole in 10^6
>>> sq km
>>>
>>> print(Pole_area)
>>>
>>> Here are the printVarSummary results:
>>>
>>> Variable: areaa
>>> Type: float
>>> Total Size: 131072 bytes
>>> 32768 values
>>> Number of Dimensions: 2
>>> Dimensions and sizes: [lat | 128] x [lon | 256]
>>> Coordinates:
>>> lat: [-88.92773..88.92773]
>>> lon: [ 0..358.5938]
>>> Number Of Attributes: 8
>>> long_name : area of each grid cell
>>> units : km^2
>>> area_total : 5.101007e+08
>>> area_lat : <ARRAY of 128 elements>
>>> rearth : 6371.22
>>> area_sphere : 5.100997e+08
>>> area_ratio : 1.000002
>>> typeConversion_op_ncl : double converted to float
>>> (0)
>>> (0) area of each grid cell: min=451.99 max=24355.4
>>>
>>>
>>> Variable: S_Area
>>> Type: double
>>> Total Size: 8 bytes
>>> 1 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [1]
>>> Coordinates:
>>> (0) 510099745.7121028
>>>
>>>
>>> Variable: Area_total
>>> Type: float
>>> Total Size: 4 bytes
>>> 1 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [1]
>>> Coordinates:
>>> (0) 5.101245e+08
>>>
>>>
>>> Variable: Pole_area
>>> Type: double
>>> Total Size: 8 bytes
>>> 1 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [1]
>>> Coordinates:
>>> (0) -0.04953257579445839
>>>
>>>
>>> Thanks,
>>> Soumik
>>>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Oct 11 13:13:48 2013