Re: Data Crossing the 180 Degree Longitude Line

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Thu, 25 Sep 2008 10:28:37 -0600

This was resolved off-line.

The final result:
====
I have attached a script with a few extra things + I commented out the

 ; find alternate max/min longitude points, if necessary

section and used the "where" statement
 http://www.ncl.ucar.edu/Document/Functions/Built-in/where.shtml

 lon2d = where(lon2d.lt.0, lon2d+360, lon2d)
===

Tony.Conrad wrote:
> 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
>

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Sep 25 2008 - 10:28:37 MDT

This archive was generated by hypermail 2.2.0 : Mon Sep 29 2008 - 13:35:43 MDT