# Re: Calculating surface area of a T42 grid

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Thu, 13 Aug 2009 16:12:03 -0600

Ummm, nobody stepped up :-(
=============================
Try the attached

Good luck
D

Gina Henderson wrote:
> Hi there,
>
> has anyone used ncl to calculate the surface area (km^2) of CCSM grid
> squares? I am looking to calculate the area of land occupied by
> certain grid squares. I have the lat and lon arrays of what I assume
> are the grid vertices and was going to use teh gc_qarea function to
> calc the area on a unit sphere and then go from there. Is there a
> better way to do this? Also, I want these grid cell areas to
> correspond to data fields that I plotted using the rasterfill cn
> funtion (on the same 64x128 grid), so does this mean I should offset
> the vertices to build the area from the 'center' of a grid cell rather
> than the corners?
>
> Any suggestions or comments would be great.
> Thanks, Gina.
>

; sfc area of spherical earth: 510,072,000 km²
; land 148,940,000 km² land (29.2 %)
; water 361,132,000 km² water (70.8 %)
;==============================================

begin
diri = "/project/cas/shea/netCDF_Misc/"
fili = "atmos.nc"

lat = f->lat
lon = f->lon
gwt = f->gw

nlat = dimsizes(lat)
mlon = dimsizes(lon)

re = 6371.220 ; [km] average radius of earth
pi = 4.0*atan(1.0)
ASPHERE= 4*pi*re^2 ; km^2 [theoretical sfc sphere]
ALAND = 148940e3 ; km^2 land area
AWATER = 361132e3 ; km^2 water area

print(" ")
print("=============> WIKIPEDIA <===========")
print("AEARTH="+ASPHERE)
print("AWATER="+AWATER)
print("ALAND ="+ALAND)

dxeq = (lon(2)-lon(1))*rr ; dx=dlon at equator [m]

dx = dxeq*cos(lat*rad) ; dx[*] at each latitude

dy = new (nlat,typeof(lat),"No_FillValue")
dy(0) = (90-abs((lat(1)+lat(0))*0.5))*rr
dy(nlat-1)= dy(0)
dy(1:nlat-2) = abs(lat(2:nlat-1)-lat(1:nlat-2))*rr

carea = dx*dy ; [*] cell area function of latitude only
CAREA = conform_dims( (/nlat,mlon/), carea, 0)

AREA = sum(CAREA) ; km^2

print(" ")
print("=============> SPHERE vs T42 <===========")
print("ASPHERE="+ASPHERE+" AREA="+AREA)
print(" ")

;+++++++++++++++++++++++++++++++++++++++++++
CAREA@_FillValue = 1e20

oro = f->ORO(0,:,:) ; "ocean (0), land (1), sea ice (2) flag"

CAREA_LAND = where(oro.eq.1, CAREA, CAREA@_FillValue)
CAREA_WATER = where(oro.eq.0, CAREA, CAREA@_FillValue)

CLAND = sum(CAREA_LAND)
CWATER = sum(CAREA_WATER)
print("=============> T42 <===========")
print("CLAND="+CLAND+" CWATER="+CWATER)
print(" ")
landDiff = ((CLAND -ALAND )/ALAND )*100 ; %
print("landDiff[%]="+ landDiff)
print(" ")
waterDiff = ((CWATER-AWATER)/AWATER)*100 ; %
print("waterDiff[%]="+ waterDiff)
end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Aug 13 2009 - 16:12:03 MDT

This archive was generated by hypermail 2.2.0 : Fri Aug 14 2009 - 23:36:52 MDT