[ncl-talk] clean edges along a satellite orbit.
Karin Meier-Fleischer
meier-fleischer at dkrz.de
Thu Jun 9 02:17:08 MDT 2016
Hi Laura,
to avoid the edge effects you can plot your data points using polymarker
(filled circle or square) or you can regrid the data using ESMF to a
rectilinear grid with a high resolution.
I've attached two example scripts reading and plotting EUMETSAT data
which seems to have the same structure like your data.
Bye,
Karin
Am 08.06.16 um 19:01 schrieb Laura Fowler:
> Hi:
>
> I am not sure if I correctly submitted the e-mail below to ncl-talk
> yesterday. So I am resubmitting it today just in case. I would really
> appreciate your help. Thanks.
> Laura
>
>
>
> I am trying to plot orbital data from the Terra Satellite, namely the
> TOA outgoing long wave radiation. The netCDF file contains 3
> one-dimensional fields that I am trying to use to plot the orbital
> data with a cylindrical-equidistant projection.
>
> LonFov(nTime): Longitude of the pixel.
> LatFov(nTime): Latitude of the pixel.
> TOAlw(nTime): Value of the outgoing long-wave radiation located at
> LonFov and LatFov.
>
> nTime is the length of my data set, so that I can plot a portion of
> the satellite orbit. To plot the data, I use
>
> res=True
> res at gsnDraw = True
> res at gsnFrame = True
> res at gsnMaximize = True
> res at gsnSpreadColors = True
>
> res at mpProjection = "CylindricalEquidistant"
> res at mpDataBaseVersion = "MediumRes"
> res at mpCenterLatF = 0.
> res at mpCenterLonF = 180.
>
> res at cnFillMode = "RasterFill"
> res at cnFillOn = True
> res at cnLinesOn = False
>
> res at cnInfoLabelOn = False
>
> res at sfXArray = LonFOV
> res at sfYArray = LatFOV
>
> wks = gsn_open_wks("pdf","test")
> plot = gsn_csm_contour_map(wks,TOAlw,res)
>
> My ncl script seems to work until the orbit widens and the edges of
> the orbits start to show some streaks because it is trying to fill
> contours (I think that this is what is happening). Is there a way that
> I can correct this. I am attaching on of my plot so that you can see
> what I am trying to do. I looked for the scripts in the Applications
> directory but did not find how to correct this issue.
>
> Thanks for your help,
> Laura
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160609/427d21e3/attachment.html
-------------- next part --------------
;-------------------------------------------------------------------------
; Plot IASI EUMETSAT original data
;
; Data: Satellite swath 3hr data (unstructured)
;
; dimensions:
; nMeas = 162417 ;
; nInstr = 1 ;
; nSens = 1 ;
; nChannel = 5 ;
; nAction = 19 ;
; nTest = 15 ;
; StringLen = 109 ;
; variables:
; float LAT(nMeas) ;
; LAT:long_name = "LATITUDE (HIGH ACCURACY)" ;
; LAT:units = "DEGREE" ;
; LAT:codetable = 5001 ;
; LAT:_FillValue_0 = 9.96921e+36f ;
; float LON(nMeas) ;
; LON:long_name = "LONGITUDE (HIGH ACCURACY)" ;
; LON:units = "DEGREE" ;
; LON:codetable = 6001 ;
; LON:_FillValue_0 = 9.96921e+36f ;
; float BT_OBS(nMeas, nChannel) ;
; BT_OBS:long_name = "Observed brightness temperature" ;
; BT_OBS:units = "K" ;
; BT_OBS:codetable = 12063 ;
; BT_OBS:_FillValue = 9.96921e+36f ;
;
;
; Plot type: polymarkers
;
; 04.02.16 kmf
;-------------------------------------------------------------------------
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
st = get_cpu_time() ;-- for calculating the elapsed CPU time
diri = "$HOME/data/EUMETSAT/IASI/"
odir = "./"
fili = "IASI_004_200907170300-200907170600.nc"
pltname = "plot_IASI_EUMETSAT_original_data_markers"
chan = 2
;-- read the data from new netCDF file
f = addfile(diri+fili,"r")
var = f->BT_OBS(:,chan)
lat = f->LAT
lon = f->LON
;-- open a workstation
wks = gsn_open_wks("png",odir+"/"+pltname)
;-- assign and retrieve color map
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
cmap = gsn_retrieve_colormap(wks)
;-- define value levels
levels = ispan(220,300,2)*1.
nlevels = dimsizes(levels)
;-- define color array
colors = span_color_indexes(cmap(3:,:),dimsizes(levels)+1) + 3
;-- set map resources
mpres = True
mpres at gsnMaximize = True ;-- maximize size of plot in window
mpres at gsnDraw = False ;-- turn off draw
mpres at gsnFrame = False ;-- turn off page advance
mpres at mpMinLatF = -90.
mpres at mpMaxLatF = 90.
mpres at mpMinLonF = -180.
mpres at mpMaxLonF = 180.
mpres at mpDataBaseVersion = "MediumRes" ;-- better map resolution
mpres at mpLandFillColor = "tan"
mpres at tiMainString = "IASI EUMETSAT: 3hr test data"
mpres at pmTickMarkDisplayMode = "Always" ;-- nicer map tickmarks
;-- create the map
map = gsn_csm_map(wks,mpres)
;-- group the variable values according to which range they fall
;-- in, and attach them to the map as a colored marker
mkres = True
mkres at gsMarkerIndex = 16 ;-- filled dot
mkres at gsMarkerSizeF = 0.001
markerid = new(nlevels+1,graphic)
do i=0,nlevels
if(i.eq.0) then ;-- first level
ii := ind(var.lt.levels(0))
else if(i.eq.nlevels) then ;-- middle levels
ii := ind(var.ge.levels(nlevels-1))
else ;-- last level
ii := ind(var.ge.levels(i-1).and.var.lt.levels(i))
end if
end if
if(.not.any(ismissing(ii))) then
mkres at gsMarkerColor = colors(i)
markerid(i) = gsn_add_polymarker(wks,map,lon(ii),lat(ii),mkres)
end if
end do
;-- draw map and the attached markers
draw(map)
frame(wks)
;-- calculate and print elapsed wall clock time
et = get_cpu_time()
print(" --> Used CPU time: "+ (et-st) + " seconds")
end
-------------- next part --------------
;-----------------------------------------------------------
; remap_IASI_EUMETSAT.ncl:
;
; Remap the IASI EUMETSAT data using ESMF to T63 gaussian grid.
;
; dimensions:
; nMeas = 162417 ;
; nInstr = 1 ;
; nSens = 1 ;
; nChannel = 5 ;
; nAction = 19 ;
; nTest = 15 ;
; StringLen = 109 ;
; variables:
; float LAT(nMeas) ;
; LAT:long_name = "LATITUDE (HIGH ACCURACY)" ;
; LAT:units = "DEGREE" ;
; LAT:codetable = 5001 ;
; LAT:_FillValue_0 = 9.96921e+36f ;
; float LON(nMeas) ;
; LON:long_name = "LONGITUDE (HIGH ACCURACY)" ;
; LON:units = "DEGREE" ;
; LON:codetable = 6001 ;
; LON:_FillValue_0 = 9.96921e+36f ;
; float BT_OBS(nMeas, nChannel) ;
; BT_OBS:long_name = "Observed brightness temperature" ;
; BT_OBS:units = "K" ;
; BT_OBS:codetable = 12063 ;
; BT_OBS:_FillValue = 9.96921e+36f ;
;
; 05.02.16 kmf
;-----------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
print("")
st = get_cpu_time() ;-- for calculating the elapsed CPU time
diri = "$HOME/data/EUMETSAT/IASI/"
fname = "IASI_004_200907170300-200907170600.nc"
odir = "./"
ofile = str_get_field(fname,1,".")+"_remap_to_T63.nc"
;-- output file: T63 grid definition
gridx = 192 ;-- number of longitude values
gridy = 96 ;-- number of latitude values
gridxinc = 1.875 ;-- longitude increment
;-- gaussian latitude values
latY = (/ -88.57217 , -86.72253 , -84.86197 , -82.99894 , \
-81.13498 , -79.27056 , -77.40589 , -75.54106 , -73.67613 , \
-71.81113 , -69.94608 , -68.08099 , -66.21587 , -64.35073 , -62.48557 , \
-60.6204 , -58.75521 , -56.89001 , -55.02481 , -53.1596 , \
-51.29438 , -49.42915 , -47.56393 , -45.69869 , -43.83346 , \
-41.96822 , -40.10298 , -38.23774 , -36.37249 , -34.50724 , \
-32.64199 , -30.77674 , -28.91149 , -27.04624 , -25.18099 , \
-23.31573 , -21.45048 , -19.58522 , -17.71996 , -15.8547 , \
-13.98945 , -12.12419 , -10.25893 , -8.393669, -6.528409, \
-4.66315 , -2.79789 , -0.9326299, 0.9326299, 2.79789 , \
4.66315 , 6.528409, 8.393669, 10.25893 , 12.12419 , \
13.98945 , 15.8547 , 17.71996 , 19.58522 , 21.45048 , \
23.31573 , 25.18099 , 27.04624 , 28.91149 , \
30.77674 , 32.64199 , 34.50724 , 36.37249 , 38.23774 , \
40.10298 , 41.96822 , 43.83346 , 45.69869 , 47.56393 , \
49.42915 , 51.29438 , 53.1596 , 55.02481 , 56.89001 , \
58.75521 , 60.6204 , 62.48557 , 64.35073 , 66.21587 , \
68.08099 , 69.94608 , 71.81113 , 73.67613 , 75.54106 , \
77.40589 , 79.27056 , 81.13498 , \
82.99894 , 84.86197 , 86.72253 , 88.57217 /)
;-- saver way to define lonX (lonX: 192 grid points, starting at -180.0 equally spaced 1.875)
lonX = ispan(0,(gridx-1),1)*gridxinc
lonX = where(lonX .lt. 180.,lonX,lonX-360.)
qsort(lonX) ;-- sort values to get -180.0-178.125
;-- create coordinate and named coordinate dimensions
latY!0 = "lat"
latY&lat = latY
latY at units = "degrees_north"
lonX!0 = "lon"
lonX&lon = lonX
lonX at units = "degrees_east"
;-- open file
if (isfilepresent(diri+"/"+fname)) then
f = addfile(diri+"/"+fname,"r")
else
print("++++ File not found !")
end if
;-- assign data (LON, LAT and BT_OBS) from source netCDF file
bt1 = f->BT_OBS(:,2)
lat1d = f->LAT(:)
lon1d = f->LON(:)
if(any(ismissing(bt1))) then
print("+++ data contains missing values ...")
end if
;-- print date to logfile to retrieve run time
date = systemfunc("date")
print("--> ESMF regridding start - "+fname+" "+date)
;-- ESMF resource settings
reopt = True ;-- ESMF resource settings
reopt at ForceOverwrite = True
reopt at WgtFileName = odir+"bt_to_T63.nc" ;-- weight file
reopt at NoPETLog = True ;-- don't generate the "PET0.RegridWeightGen.log" file
reopt at RemoveSrcFile = True ;-- remove Src SCRIP grid destination files
reopt at RemoveDstFile = True ;-- remove Dst SCRIP grid destination files
reopt at SrcFileName = odir+"bt2d_ESMF.nc" ;-- output file
reopt at DstTitle = "Swath T63 resolution" ;-- destination file title
reopt at DstGridLat = latY ;-- destination lat grid values
reopt at DstGridLon = lonX ;-- destination lon grid values
reopt at DstRegional = True
; reopt at Debug = True
reopt at SrcRegional = True
reopt at DstFileName = "Swath_T63_SCRIP.nc" ;-- output SCRIP file
;-- open a workstation
wks = gsn_open_wks("png",odir+"plot_IASI_regridded_contour_bilinear_large")
res = True
res at gsnAddCyclic = False
res at gsnMaximize = True
res at cnFillOn = True
res at cnFillMode = "RasterFill"
res at cnLinesOn = False
res at cnInfoLabelOn = False
res at cnLineLabelsOn = False
res at cnLevelSelectionMode = "ManualLevels"
res at cnMinLevelValF = 220.
res at cnMaxLevelValF = 300.
res at cnLevelSpacingF = 2.
res at cnFillPalette = "BlAqGrYeOrReVi200"
res at gsnFrame = False
res at gsnDraw = False
cnres = res
cnres at gsnRightString = ""
cnres at gsnLeftString = ""
cnres at lbLabelBarOn = False
mpres = res
mpres at mpOutlineOn = True
mpres at mpGridAndLimbOn = True
mpres at mpGridLineColor = "grey30"
mpres at mpGridLineThicknessF = 1
mpres at mpGridLineDashPattern = 0
mpres at mpGridAndLimbDrawOrder = "PostDraw"
mpres at tiMainString = "IASI EUMETSAT - ESMF regridded to T63"
;-- regrid the data in rectangular blocks
lats = ispan(-90,90,45) ;-- result rough but the best
lons = ispan(-180,180,90) ;-- result rough but the best
nlats = dimsizes(lats)
nlons = dimsizes(lons)
first_segment = True
do nlt=0,nlats-2
do nln=0,nlons-2
minlat = lats(nlt)
maxlat = lats(nlt+1)
minlon = lons(nln)
maxlon = lons(nln+1)
;-- find all the lat/lon points in this rectangular block
ii := ind(lat1d.ge.minlat-1.and.lat1d.le.maxlat+1.and.lon1d.ge.minlon-1.and.lon1d.le.maxlon+1)
print("Range lat="+minlat+":"+maxlat+", lon="+minlon+":"+maxlon)
if(any(ismissing(ii)).or.dimsizes(ii).le.3) then
print("No lat/lon values found in this range.")
continue
else
print(dimsizes(ii) + " values found in this range.")
end if
;-- ESMF regridding
reopt at SrcGridLat := lat1d(ii) ;-- source lat grid values
reopt at SrcGridLon := lon1d(ii) ;-- source lon grid values
bt1_regrid = ESMF_regrid(bt1(ii), reopt) ;-- take some time
;---Create plots of each subset of data.
if(first_segment) then
plot_regrid = gsn_csm_contour_map(wks,bt1_regrid,mpres)
first_segment = False
else
overlay_regrid = gsn_csm_contour(wks,bt1_regrid,cnres)
overlay(plot_regrid,overlay_regrid)
end if
end do
end do
draw(plot_regrid)
frame(wks)
;-- print date to logfile to retrieve run time
et = get_cpu_time() ;-- for calculating the elapsed CPU time
cputime = et-st
print("----> CPU time : "+cputime+" s")
end
More information about the ncl-talk
mailing list