| co-kringing [message #75662] |
Thu, 07 April 2011 06:49 |
ggggg
Messages: 2
Registered: January 2011
|
Junior Member |
|
|
hallo all,
I tried to do some kriging with irregularly sampled precipation
measurments, worked not bad (have a look at the code) know i want do
do some cokriging with precipation and height. is that possible in
idl? which procedure does help me? couldn't find anything!
as input i have 4 one-dimensional arrays with
x = stationlist.lon ;longitudes fltarr(600)
y = stationlist.lat ;latitudes fltarr(600)
data2= total(data,2) ;precipation data fltarr(600)
data3 ;height data fltarr(600)
hope someone can help me
thanks
M.S.
####code##########
pro example7griddata
restore,filename='h:\schule\bachelorarbeit\daten3km.sav'
maps
;für alle modelldaten########################
x = stationlist.lon ;longitudes fltarr(600)
y = stationlist.lat ;latitudes fltarr(600)
data2= total(data,2) ;precipation data fltarr(600)
data3 ;height data fltarr(600)
;###########################3
; Grid the Data and Display the Results:
; Preprocess and sort the data. GRID_INPUT will
; remove any duplicate locations.
GRID_INPUT, x, y, data2, xSorted, ySorted, dataSorted
; Initialize the grid parameters.
gridSize = [100, 100]
minxy=[5.9,45.5] ;minxy=[minx,miny] boundaries of the contour plot,
bottom left corner
maxxy=[10.6,48];maxxy=[maxx,maxy] upper right corner
; Use the equation of a straight line and the grid parameters to
; determine the x of the resulting grid.
slope = (MAXxy[0] - MINxy[0])/(gridSize[0] - 1)
intercept = Minxy[0]
xGrid = (slope*FINDGEN(gridSize[0])) + intercept
; Use the equation of a straight line and the grid parameters to
; determine the y of the resulting grid.
slope = (MAXxy[1] - MINxy[1])/(gridSize[1] - 1)
intercept = MINxy[1]
yGrid = (slope*FINDGEN(gridSize[1])) + intercept
DEVICE, DECOMPOSED = 0
loadct,1
; Grid the data with faulting using the Radial Basis Function
; method.
grid = GRIDDATA(xSorted, ySorted, dataSorted, $
DIMENSION = gridSize, METHOD = 'kriging', $
MISSING = MIN(dataSorted))
col=findgen(255)
col=reverse(col)
; Display grid results.
CONTOUR, BYTSCL(grid), xGrid, YGrid,/overplot, /follow,/XSTYLE, /
YSTYLE, nlevels=n_elements(col),C_COLORS =col,/fill
MAP_CONTINENTS, /COUNTRIES, COLOR=white, MLINETHICK=2, /hires
MAP_CONTINENTS, /COASTS, COLOR=255, /hires
end
Pro maps
maplimit = [45.7, 5.9, 47.8, 10.6] ;NMM22 europe
; Handle TrueColor displays:
DEVICE, DECOMPOSED=0
; Load discrete color table:
tek_color
; Match color indices to colors we want to use:
black=0 & white=1 & red=2
green=3 & dk_blue=4 & lt_blue=5
; Set up an orthographic projection centered over the north
; Atlantic.Fill the hemisphere with dark blue. Specify black
; gridlines:
;MAP_SET, /ORTHO, 40, -30, 23, /ISOTROPIC, $
; /HORIZON, E_HORIZON={FILL:1, COLOR:dk_blue}, $
; /GRID, COLOR=black
MAP_SET, /ORTHOGRAPHIC, 48, 7, limit = maplimit,/iso, $
/GRID, LONDEL=15, LATDEL=15,color=1, title ='stationen',/hires
; Fill the continent boundaries with solid white:
MAP_CONTINENTS, /FILL_CONTINENTS, COLOR=1,/hires
; Overplot coastline data:
MAP_CONTINENTS, /COASTS, COLOR=0, /hires
; Add rivers, in light blue:
;MAP_CONTINENTS, /RIVERS, COLOR=lt_blue
; Show national borders:
MAP_CONTINENTS, /COUNTRIES, COLOR=0, MLINETHICK=2, /hires
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
