Re: Transform data to a different map projection

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri, 6 Oct 2006 19:20:22 -0600 (MDT)

> I want to transform a data of Lambert Conformal Format to Cylindrical Equidistant projection (regular lat/lon grid). More inflormation on this format can be found here: http://www.cdc.noaa.gov/NARR/format.html. Is there a simple way to transform it in NCL? Can I plot it in Cylindrical Equidistant projection without transforming?
>
> The description of the Lambert Conformal Format data is:
>
> dimensions:
> time = UNLIMITED ; // (365 currently)
> y = 277 ;
> x = 349 ;
> nbnds = 2 ;
> variables:
> float lat(y, x) ;
> lat:long_name = "latitude coordinate" ;
> lat:units = "degrees_north" ;
> lat:axis = "Y" ;
> lat:coordinate_defines = "point" ;
> lat:standard_name = "latitude" ;
> float lon(y, x) ;
> lon:units = "degrees_east" ;
> lon:long_name = "longitude coordinate" ;
> lon:axis = "X" ;
> lon:coordinate_defines = "point" ;
> lon:standard_name = "longitude" ;
> float x(x) ;
> x:long_name = "eastward distance from southwest corner of domain in projection coordinates" ;
> x:units = "m" ;
> x:standard_name = "projection_x_coordinate" ;
> float y(y) ;
> y:long_name = "northward distance from southwest corner of domain in projection coordinates" ;
> y:units = "m" ;
> y:standard_name = "projection_y_coordinate" ;
> int Lambert_Conformal ;
> Lambert_Conformal:grid_mapping_name = "lambert_conformal_conic" ;
> Lambert_Conformal:standard_parallel = 50., 50. ;
> Lambert_Conformal:longitude_of_central_meridian = -107. ;
> Lambert_Conformal:latitude_of_projection_origin = 50. ;
> Lambert_Conformal:false_easting = 5632642.22547 ;
> Lambert_Conformal:false_northing = 4612545.65137 ;
> double time(time) ;
> time:units = "hours since 1800-1-1 00:00:0.0" ;
> time:long_name = "analysis time" ;
> time:axis = "T" ;
> time:standard_name = "time" ;
> time:coordinate_defines = "start" ;
> time:delta_t = "0000-00-01 00:00:00" ;
> time:actual_range = 1569072., 1577808. ;
> time:avg_period = "0000-00-01 00:00:00" ;
> double time_bnds(time, nbnds) ;
> time_bnds:long_name = "Time Boundaries" ;
> short acpcp(time, y, x) ;
> acpcp:units = "kg/m^2" ;
> acpcp:long_name = "Daily accumulated convective precipitation at Surface" ;
> acpcp:unpacked_valid_range = 0.f, 60.f ;
> acpcp:precision = 3s ;
> acpcp:actual_range = 0.000995636f, 32.766f ;
> acpcp:missing_value = 32766s ;
> acpcp:add_offset = 32.766f ;
> acpcp:scale_factor = 0.001f ;
> acpcp:valid_range = -32766s, 27235s ;
> acpcp:_FillValue = -32767s ;
> acpcp:GRIB_name = "ACPCP" ;
> acpcp:GRIB_id = 63 ;
> acpcp:var_desc = "convective precipitation accumulation" ;
> acpcp:standard_name = "convective_precipitation_amount" ;
> acpcp:level_desc = "Surface" ;
> acpcp:dataset = "NARR Daily Averages" ;
> acpcp:statistic = "Mean" ;
> acpcp:parent_stat = "Individual Obs" ;
> acpcp:grid_mapping = "Lambert_Conformal" ;
> acpcp:coordinates = "lat lon" ;
> acpcp:cell_methods = "time: mean (of 8 3-hourly values in one day)" ;
>
> // global attributes:
> :standardpar1 = 50. ;
> :standardpar2 = 50.000001 ;
> :centerlon = -107. ;
> :centerlat = 50. ;
> :latcorners = 1.000001f, 0.897945f, 46.3544f, 46.63433f ;
> :loncorners = -145.5f, -68.32005f, -2.569891f, 148.6418f ;
> :stream = "s1" ;
>
>
>
> -----------------------------

[1] Click the "Applications" and under "Datasets" click NARR

     http://www.ncl.ucar.edu/Applications/narr.shtml

[2] If you wantto just plot the NARR on a cylindrical projection
     then the following modification of narr_2.ncl might work

     NOT TESTED

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
;*******************************************
; open file and read in data
;*******************************************
   f = addfile ("foo.nc", "r")
   lat2d = f->y
   lon2d = f->x

   print(" min="+min(lat2d)+" max="+max(lat2d))
   print(" min="+min(lon2d)+" max="+max(lon2d))

   time = f->time
   date = ut_calendar (time, -2)
   date_at_units = "yyyymmdd"
  ;print(date)

   x = f->acpcp
   printVarSummary(x)
   print(" min="+min(x)+" max="+max(x))

;********************************************
; create plots
;********************************************
   wks = gsn_open_wks ("ps", "narr") ; open workstation
   gsn_define_colormap (wks,"gui_default") ; choose color map

   res = True ; plot mods desired for original grid
   res_at_cnFillOn = True ; color fill
   res_at_cnLinesOn = False ; no contour lines
   res_at_cnLineLabelsOn = False ; no contour labels
   res_at_gsnSpreadColors = True ; use total colormap
   res_at_cnInfoLabelOn = False ; no contour info label

   res_at_mpGridAndLimbOn = True
   res_at_mpGridLineDashPattern = 2 ; lat/lon lines as dashed
   res_at_pmTickMarkDisplayMode = "Always" ; turn on tickmarks
   res_at_tmXTOn = False

   res_at_gsnAddCyclic = False ; regional data

   dimlc = dimsizes(lat2d)
   nlat = dimlc(0)
   mlon = dimlc(1)

   res_at_mpLimitMode = "Corners" ; choose range of map
   res_at_mpLeftCornerLatF = lat2d(0,0)
   res_at_mpLeftCornerLonF = lon2d(0,0)
   res_at_mpRightCornerLatF = lat2d(nlat-1,mlon-1)
   res_at_mpRightCornerLonF = lon2d(nlat-1,mlon-1)
   res_at_mpProjection = "CylindricalEquidistant"

   nt = 0
   res_at_gsnCenterString = date(nt)
   plot = gsn_csm_contour_map(wks,x(nt,:,:),res) ; Draw original grid

   end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Oct 06 2006 - 19:20:22 MDT

This archive was generated by hypermail 2.2.0 : Mon Oct 09 2006 - 09:24:50 MDT