begin
;y = cbinread("./N22E113.hgt",-1,"short")
;printVarSummary(y)
setfileoption ("bin", "ReadByteOrder", "BigEndian")
demdata = cbinread("./N22E113.hgt",(/3601,3601/),"short")
;print(demdata(:,0))
;print(max(demdata))
;print(min(demdata))
lat_s = 22.00
lat_n = 23.00
lon_w = 113.0
lon_e = 114.0
;nlat = 3601
;nlon = 3601
;lat1d = fspan(lat_s,lat_n,nlat) ;---结果出现坐标间距可变bug
;lon1d = fspan(lon_w,lon_e,nlon)
dxy = 1d0/3600 ;---双精度---
fxy = fspan(0, 3600*dxy, 3601)
lon1d = fxy + lon_w
lat1d = lat_n - fxy ;---y reverse
demdata!0 = "lat"
demdata!1 = "lon"
demdata&lat = lat1d
demdata&lon = lon1d
demdata&lat@units = "degrees_north"
demdata&lon@units = "degrees_east"
;---write---
system("/bin/rm -f PRDdem_o.nc") ; remove any pre-existing file
ncdf = addfile("PRDdem_o.nc" ,"c") ; open output netCDF file
ncdf->z = demdata
end
