Data Crossing the 180 Degree Longitude Line

Hello Again:

This script (reproduced below), which maps sea surface temperature (sst)
data on a Mercator projection, works as written. But I keep wondering
if there is a better way to do this.

The issue is that the sst data goes from roughly longitude -166 degrees
(or 166 degrees West) to roughly +166 degrees (or 166 degrees East),
crossing the 180 degree longitude line. In other words, the data goes
from -166 to -180 degrees, then from +166 to +180 degrees. Since the
maximum longitude is +180 and the minimum longitude is -180, whenever
I've tried mapping the sst points NCL will draw the entire globe with a
tiny square of data centering on the 180 degree longitude line [[the sst
data goes from latitude 4 to 27 degrees North]]. The way to stop this
is to employ two do-loops in the script (found under the comment "find
alternate max/min longitude points, if necessary"). These two do-loops
find the maximum negative longitude (-166) and the minimum positive
longitude (+166) and those are used as max/min longitudes instead of
-180/+180. But--is there a better way to do this?

---Thanx
---Tony Conrad
---NSOF/NOAA
---Suitland, Maryland, USA

;*********************************
; hawaiim_4.ncl
; Modis data on a Mercator projection.
;*********************************
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
;*********************************
; read in HDF file
;*********************************
  f1 =
addfile("/usr/data/modisRedo/223-Day/A2008223012500.L2_LAC_SST.hdf", "r")
;*************************************************************
; read in sst data
;*************************************************************
  ch1_int = f1->sst
  ch1 = ch1_int*ch1_int_at_slope(0) ; orig is integer, make float

;*********************************
; read in second HDF file
;*********************************
  f2 =
addfile("/usr/data/modisRedo/223-Day/MYD03.A2008223.0125.005.2008223165935.hdf",
"r")
;*************************************************************
; read in navigation data about the sst data
;*************************************************************
  lat2d = f2->Latitude
  lon2d = f2->Longitude

;*********************************
; read in third HDF file
;*********************************
  f3 =
addfile("/usr/data/modisRedo/223-Day/MYD35_L2.A2008223.0125.005.2008223174653.hdf",
"r")
;*************************************************************
; read in cloud mask data about the sst data
;*************************************************************
  cloud = f3->Cloud_Mask(0,:,:)

;*************************************************************
; isolate least significant bits (cloud mask) in cloud mask data
;*************************************************************
  c_bits = dim_gbits( cloud, 5, 3, 5, 1354 ) ; isolate three bits
;*************************************************************
; isolate most significant bits (water/land bits) in cloud mask data
;*************************************************************
  w_bits = dim_gbits( cloud, 0, 2, 6, 1354 ) ; isolate two bits

;*************************************************************
; mask out cloudy sst data and sst data over land
;*************************************************************
  ch2 = mask( ch1, (c_bits.ne.1 .and. w_bits.eq.0), True )

  ch2_at_lat2d = lat2d ; assign nav data to sst data
  ch2_at_lon2d = lon2d
;*************************************************************
; assign named dimensions
;*************************************************************
  ch2_at_long_name = "Sea Surface Temp"
  ch2_at_units = "degrees Celsius"

;*************************************************************
; find alternate max/min longitude points, if necessary
;*************************************************************
  if(( max(lon2d) .gt. 179.9 ) .and. ( min(lon2d) .lt. -179.9 )) then
     lon1d = ndtooned( lon2d ) ; change lon array from 2d to 1d
     lon1d_size = dimsizes( lon1d ) ; get size of lon array
     maxlon = 180.0 ; start at max lon (max East)
     do q = 0, lon1d_size-1 ; step thru lon array
        if( lon1d(q) .lt. 0.0 ) then ; toss out all negative values
           continue
        end if
        if( lon1d(q) .lt. maxlon ) then ; continue replacing max with
           maxlon = lon1d(q) ; smaller value--should end up
with
        end if ; smallest positive longitude
angle
     end do
     minlon = -180.0 ; start at min lon (max West)
     do q = 0, lon1d_size-1 ; step thru lon array
        if( lon1d(q) .gt. 0.0 ) then ; toss out all positive values
           continue
        end if
        if( lon1d(q) .gt. minlon ) then ; continue replacing min with
           minlon = lon1d(q) ; larger value--should end up with
        end if ; largest negative longitude angle
     end do
  end if

;*************************************************************
; create plot
;*************************************************************
  wks = gsn_open_wks("x11", "hawaiim")
  gsn_define_colormap(wks, "rainbow+white+gray") ; choose colormap
 
  res = True
  res_at_cnFillOn = True ; color Fill
  res_at_cnFillMode = "CellFill" ; fill mode
  res_at_cnLinesOn = False ; Turn off contour lines
  res_at_cnLineLabelsOn = False ; Turn off contour lines
 
  res_at_cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  res_at_cnMinLevelValF = -50.0 ; set min contour level
  res_at_cnMaxLevelValF = 30.0 ; set max contour level
  res_at_cnLevelSpacingF = 1.0 ; set contour spacing
 
  res_at_gsnSpreadColors = True ; use full colormap
  res_at_gsnSpreadColorStart = 42 ; start at color 42
  res_at_gsnSpreadColorEnd = 237 ; end at color 237
  res_at_lbLabelStride = 10 ; every 10nth label bar label
  res_at_lbOrientation = "Vertical" ; vertical label bar
 
  res_at_mpProjection = "Mercator" ; choose projection
  res_at_mpLimitMode = "Corners" ; method to zoom
  res_at_mpLeftCornerLatF = min(lat2d) ; choose lat/lon map limits
  res_at_mpRightCornerLatF = max(lat2d)
  if(( max(lon2d) .gt. 179.9 ) .and. ( min(lon2d) .lt. -179.9 )) then
     res_at_mpLeftCornerLonF = maxlon
     res_at_mpRightCornerLonF = minlon
  else
     res_at_mpLeftCornerLonF = min(lon2d)
     res_at_mpRightCornerLonF = max(lon2d)
  end if
  res_at_mpCenterLonF = 180 ; change map center
  res_at_gsnMaximize = True ; maximize (without
                                            ; overflowing) the plot area

;*************************************************************
; turn on lat/lon labeling
;*************************************************************
  res_at_pmTickMarkDisplayMode = "Always" ; turn on tickmarks
  res_at_mpGridAndLimbOn = True ; turn on lat/lon lines
  res_at_mpGridLatSpacingF = 5 ; spacing for lat lines
  res_at_mpGridLonSpacingF = 10 ; spacing for lon lines

  res_at_mpFillDrawOrder = "PreDraw" ; draw land last
  res_at_mpLandFillColor = 191 ; color for land
  res_at_mpOceanFillColor = 32 ; color for ocean
  res_at_cnMissingValFillColor = -1 ; fill color = transparent
  res_at_mpDataBaseVersion = "Ncarg4_1" ; higher res coastline
  res_at_gsnAddCyclic = False ; needed b/c data not global
 
  res_at_tiMainString = "Modis data on a Mercator projection"
 
  plot = gsn_csm_contour_map(wks, ch2, res) ; create plot
 
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Sep 23 2008 - 11:29:43 MDT