[ncl-talk] Topography with GTOPO30 DATA

[Rick Brownrigg]

Hi Gerardo,

Not sure which hemisphere you are going for, but your data is from the
eastern hemisphere (60 .. 100)  and your map specification is from the
western  (-100 ... -60)  The elevation data in you plot does indeed
resemble the south china sea area, whereas the country outlines are clearly
parts of south, central, and north america.

hope that helps...
Rick


On Mon, Nov 7, 2016 at 12:00 PM, Gerardo Montoya <gemonga at gmail.com> wrote:

> Dear users, I introduced some changes in the script. topo_1.ncl, in order
> to run it for high resolution elevation data (GTOPO30, tile E100N40.DEM).
> However, some things I'm doing wrong or I'm missing, because the displayed
> topography is not so good (see attached image). Also, the maximum height in
> the domain appears to be 7213, i.e., about 273m higher than the maximum
> peak in the Americas ( 6960m the Aconcagua´s peak in Argentina). I´m
> using NCL version 6.3.0 and I´m running this program in a LENOVO Z580
> laptop with a i7 3520M processor.
>
>
>
> I guess that the problem may be in the type of the variable “elev”,
> however still now I don´t know how to solve it.  Any help for correcting
> the script is greatly appreciated. The script is:
>
>
> ;================================================
>
> ; 2016-11-7  topo1km2V_0.ncl. It is a modified version of topo_1.ncl
>
> ;================================================
>
> undef("read_elev_data")
>
> function read_elev_data(topo_file)
>
> local nlat, nlon, topo_file, lat, lon
>
> begin
>
> ;---Read data as a straight binary file
>
>   nlat =6000
>
>   nlon = 4800
>
>   setfileoption("bin","ReadByteOrder","BigEndian")
>
>   elev = cbinread(topo_file,(/nlat,nlon/),"short")
>
> ;*******begin changes****************************************
>
> ;---Create 1D coordinate arrays
>
> ;  lat       = fspan(90,-90,nlat)
>
> ;  lon       = fspan(0,360,nlon)
>
> res=0.00833333333333            ; this is the resolution of the GTOPO30
> data.
>
> ulat=39.99583333333333
>
> ulon=-99.99583333333333
>
> lat = ulat - ispan(0,(nlat-1),1)*res    ;from 39.99583333333333 to
> -9.99583333331
>
> lon = ulon+ispan(0,(nlon-1),1)*res   ;from -99.99583333333333 to
> -60.00416666668
>
>   lat at _FillValue=999.99
>
>   lon at _FillValue=999.99
>
> ;*****end changues****************************************************
> *************
>
>   lat!0     = "lat"
>
>   lon!0     = "lon"
>
>   lat at units = "degrees_north"
>
>   lon at units = "degrees_east"
>
>   lat&lat   = lat
>
>   lon&lon   = lon
>
>
> ;---Attach the coordinate arrays
>
>   elev!0    = "lat"
>
>   elev!1    = "lon"
>
>   elev&lat  = lat
>
>   elev&lon  = lon
>
>
>
>   return(elev)
>
> end
>
>
> ;----------------------------------------------------------------------
>
> ; Main code
>
> ;----------------------------------------------------------------------
>
> begin
>
>   wks = gsn_open_wks("png","topo1km2V_0")   ; send graphics to PNG file
>
> ;**************begin changues**************************
>
> setvalues NhlGetWorkspaceObjectId()
>
>  "wsMaximumSize" : 300000000
>
>  end setvalues
>
> ; the tile, E100N40.DEM was downloaded from:
>
> ; https://dds.cr.usgs.gov/srtm/version2_1/SRTM30/e100n40/
>
>   elev = read_elev_data("/usr/local/Topog/E100N40.DEM")
>
>   printVarSummary(elev)
>
>   printMinMax (elev, False)
>
> ;*************************end changues*****************
>
> ;---Set some plot options
>
>   res                    = True
>
>   res at gsnMaximize        = True             ; maximize plot in frame
>
>   res at cnFillOn           = True             ; turn on contour fill
>
>   res at cnLevelSpacingF    = 125              ; NCL picks 2000
>
>   res at cnFillMode         = "RasterFill"     ; much faster than AreaFill
>
>   res at cnLinesOn          = False            ; turn off contour lines
>
>   res at cnLineLabelsOn     = False            ; turn off line labels
>
>   res at cnInfoLabelOn      = False            ; turn off info label
>
>   res at lbBoxLinesOn       = False            ; turn off labelbar box lines
>
>   res at gsnAddCyclic       = False            ; don't add longitude cyclic
> point
>
>   res at mpFillOn           = False            ; turn off map fill
>
>   res at tiMainString       = "topo1km2V_0"     ; main title
>
>   res at pmLabelBarWidthF   = 0.8              ; default is too short
>
> ;*************begin changues*************************
>
> ; set plot domain for the tile, E100N40.DEM
>
> res at mpMinLatF           = -10.0
>
> res at mpMaxLatF           =  40.0
>
> res at mpMinLonF           = -100.0
>
> res at mpMaxLonF           = -60.0
>
> ;************end changues***************************
>
>   plot = gsn_csm_contour_map(wks,elev,res)
>
> end
>
> Thanks,
>
> Gerardo Montoya
>
> full professor (retired), Universidad Nacional de Colombia
>
> _______________________________________________
> 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/20161107/fe05c385/attachment.html