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 %)
;==============================================

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
  diri = "/project/cas/shea/netCDF_Misc/"
  fili = "atmos.nc"
  f = addfile (diri+fili, "r")

  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)
  rad = pi/180.0
  rr = re*rad
  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
  areaDiff = ((AREA-ASPHERE)/ASPHERE)*100 ; %

  print(" ")
  print("=============> SPHERE vs T42 <===========")
  print("ASPHERE="+ASPHERE+" AREA="+AREA)
  print("areaDiff[%]="+ areaDiff)
  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