[ncl-talk] Topography with GTOPO30 DATA

[Gerardo Montoya]

Thanks Rick. it was the error. I did a silly mistake downloading the tile  "
E100N40.DEM" instead of "W100N40.DEM". Now, with downloaded W100N40.DEM tile,
the graphic looks much better. Thanks newly to you and  also to the NCL
team for your fruitful work.

2016-11-07 18:07 GMT-03:00 Rick Brownrigg <brownrig at ucar.edu>:

> 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/cb3595d2/attachment.html