diff --git a/Data/NetCDF.hs b/Data/NetCDF.hs --- a/Data/NetCDF.hs +++ b/Data/NetCDF.hs @@ -46,7 +46,7 @@ , module Data.NetCDF.Metadata , NcStorable (..) , IOMode (..) - , openFile, createFile, closeFile + , openFile, createFile, syncFile, closeFile , withReadFile, withCreateFile , get1, get, getA, getS , put1, put, putA, putS @@ -98,6 +98,10 @@ chk $ nc_enddef ncid return $ NcInfo n ds vs as ncid varids +-- | Sync a NetCDF file. +syncFile :: NcInfo NcWrite -> IO () +syncFile (NcInfo _ _ _ _ ncid _) = void $ nc_sync ncid + -- | Close a NetCDF file. closeFile :: NcInfo a -> IO () closeFile (NcInfo _ _ _ _ ncid _) = void $ nc_close ncid @@ -216,7 +220,7 @@ write1Var ncid dimidmap (NcVar n t dims as) = do let dimids = map ((dimidmap M.!) . ncDimName) dims varid <- chk $ nc_def_var ncid n (fromEnum t) (length dims) dimids - forM_ (M.toList as) (write1Attr ncid varid) + forM_ (M.toList as) $ write1Attr ncid varid return varid -- | Read an attribute from a NetCDF variable with error handling. @@ -239,21 +243,23 @@ -- | Write an attribute to a NetCDF variable with error handling. writeAttr :: Int -> Int -> String -> NcAttr -> Access () -writeAttr nc var n (NcAttrByte v) = writeAttr' nc var n id v nc_put_att_uchar -writeAttr nc var n (NcAttrChar v) = writeAttr' nc var n id v nc_put_att_text +writeAttr nc var n (NcAttrByte v) = + writeAttr' nc var n NcByte id v nc_put_att_uchar +writeAttr nc var n (NcAttrChar v) = + writeAttr' nc var n NcChar id v nc_put_att_text writeAttr nc var n (NcAttrShort v) = - writeAttr' nc var n fromIntegral v nc_put_att_short + writeAttr' nc var n NcShort fromIntegral v nc_put_att_short writeAttr nc var n (NcAttrInt v) = - writeAttr' nc var n fromIntegral v nc_put_att_int + writeAttr' nc var n NcInt fromIntegral v nc_put_att_int writeAttr nc var n (NcAttrFloat v) = - writeAttr' nc var n realToFrac v nc_put_att_float + writeAttr' nc var n NcFloat realToFrac v nc_put_att_float writeAttr nc var n (NcAttrDouble v) = - writeAttr' nc var n realToFrac v nc_put_att_double + writeAttr' nc var n NcDouble realToFrac v nc_put_att_double -- | Helper function for attribute writeing. -writeAttr' :: Int -> Int -> String -> (a -> b) -> [a] +writeAttr' :: Int -> Int -> String -> NcType -> (a -> b) -> [a] -> (Int -> Int -> String -> Int -> [b] -> IO Int) -> Access () -writeAttr' nc var n conv vs wf = chk $ wf nc var n (length vs) (map conv vs) +writeAttr' nc var n t conv vs wf = chk $ wf nc var n (fromEnum t) (map conv vs) -- | Apply COARDS value scaling. diff --git a/Data/NetCDF/HMatrix.hs b/Data/NetCDF/HMatrix.hs --- a/Data/NetCDF/HMatrix.hs +++ b/Data/NetCDF/HMatrix.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} {-# OPTIONS_GHC -fno-warn-orphans #-} -- | NetCDF store instance for HMatrix vectors and matrices. module Data.NetCDF.HMatrix @@ -10,6 +10,7 @@ import Data.NetCDF.Store import qualified Numeric.Container as C +import Foreign.C import Data.Packed.Foreign import Data.Packed.Development @@ -42,3 +43,11 @@ in HColumnMajorMatrix $ matrixFromVector ColumnMajor (d `div` c) (last s) $ unsafeFromForeignPtr p 0 (Prelude.product s) smap f (HColumnMajorMatrix m) = HColumnMajorMatrix $ C.mapMatrix f m + +instance C.Element CShort +instance C.Element CInt +instance C.Element CFloat +instance C.Element CDouble + +instance C.Container C.Vector CFloat +instance C.Container C.Vector CDouble diff --git a/Data/NetCDF/Raw/Attributes.chs b/Data/NetCDF/Raw/Attributes.chs --- a/Data/NetCDF/Raw/Attributes.chs +++ b/Data/NetCDF/Raw/Attributes.chs @@ -5,13 +5,36 @@ module Data.NetCDF.Raw.Attributes where import Data.Char -import Control.Monad (liftM) import Data.NetCDF.Raw.Utils #include -nc_put_att :: (Storable a, Storable b) => +-- nc_put_att_text has to be treated seperately since it has a +-- different calling sequence to all the other +-- nc_put_att_... functions. + +-- int nc_put_att_text(int ncid, int varid, const char *name, +-- size_t len, const char *op); +nc_put_att_text :: Int -> Int -> String -> Int -> String -> IO Int +nc_put_att_text nc var name xtype v = do + let ncid = fromIntegral nc + varid = fromIntegral var + ncxtype = fromIntegral xtype + ncsize = fromIntegral $ length v + let tv = map convChar v + withCString name $ \namep -> do + withCString tv $ \vp -> do + res <- nc_put_att_text'_ ncid varid namep ncsize vp + return $ fromIntegral res + where convChar c + | isAscii c = c + | otherwise = ' ' +foreign import ccall safe "Data/NetCDF/Raw.chs.h nc_put_att_text" + nc_put_att_text'_ :: CInt -> CInt -> CString -> CULong + -> Ptr CChar -> IO CInt + +nc_put_att :: (Storable a, Storable b, Show b) => (CInt -> CInt -> CString -> CInt -> CULong -> Ptr b -> IO CInt) -> (a -> b) -> Int -> Int -> String -> Int -> [a] -> IO Int nc_put_att cfn conv nc var name xtype v = do @@ -19,21 +42,11 @@ varid = fromIntegral var ncxtype = fromIntegral xtype ncsize = fromIntegral $ length v + let tv = map conv v withCString name $ \namep -> do - (withArray . liftM conv) v $ \vp -> do + withArray tv $ \vp -> do res <- cfn ncid varid namep ncxtype ncsize vp return $ fromIntegral res - --- int nc_put_att_text(int ncid, int varid, const char *name, nc_type xtype, --- size_t len, const char *op); -nc_put_att_text :: Int -> Int -> String -> Int -> String -> IO Int -nc_put_att_text = nc_put_att nc_put_att_text'_ convChar - where convChar c - | isAscii c = fromIntegral (ord c) - | otherwise = 32 -foreign import ccall safe "Data/NetCDF/Raw.chs.h nc_put_att_text" - nc_put_att_text'_ :: CInt -> CInt -> CString -> CInt -> CULong - -> Ptr CChar -> IO CInt -- int nc_put_att_uchar(int ncid, int varid, const char *name, nc_type xtype, -- size_t len, const unsigned char *op); diff --git a/examples/example1.hs b/examples/example1.hs new file mode 100644 --- /dev/null +++ b/examples/example1.hs @@ -0,0 +1,64 @@ +module Main where + +import Prelude hiding (length, maximum, minimum, sum) +import Control.Applicative ((<$>)) +import Data.NetCDF +import Data.NetCDF.Vector +import qualified Data.Map as M +import Foreign.C +import Data.Vector.Generic hiding ((++), map) +import qualified Data.Vector.Storable as SV + +type SVRet a = IO (Either NcError (SV.Vector a)) + +main :: IO () +main = do + -- Open file and access some basic metadata. + Right nc <- openFile "z500-short.nc" + putStrLn $ "Name: " ++ ncName nc + putStrLn $ "Dims: " ++ show (M.keys $ ncDims nc) + putStr $ unlines $ map (\(n, s) -> " " ++ n ++ ": " ++ s) $ + M.toList $ flip M.map (ncDims nc) $ + \d -> show (ncDimLength d) ++ if ncDimUnlimited d then " (UNLIM)" else "" + putStrLn $ "Vars: " ++ show (M.keys $ ncVars nc) + putStrLn $ "Global attributes: " ++ show (M.keys $ ncAttrs nc) + + -- Get dimension sizes. + let Just ntime = ncDimLength <$> ncDim nc "time" + Just nlat = ncDimLength <$> ncDim nc "latitude" + Just nlon = ncDimLength <$> ncDim nc "longitude" + + -- Get longitude variable: metadata, then all values, then a slice + -- of values. + let (Just lonvar) = ncVar nc "longitude" + Right lon <- get nc lonvar :: SVRet CFloat + let mlon = mean lon + putStrLn $ "longitude: " ++ show lon ++ " -> " ++ show mlon + Right lon2 <- getS nc lonvar [0] [72] [2] :: SVRet CFloat + let mlon2 = mean lon2 + putStrLn $ "longitude (every 2): " ++ show lon2 ++ " -> " ++ show mlon2 + + -- Get latitude variable: metadata, then all values, then a slice of + -- values. + let (Just latvar) = ncVar nc "latitude" + Right lat <- get nc latvar :: SVRet CFloat + let mlat = mean lat + putStrLn $ "latitude: " ++ show lat ++ " -> " ++ show mlat + Right lat2 <- getS nc latvar [0] [37] [2] :: SVRet CFloat + let mlat2 = mean lat2 + putStrLn $ "latitude (every 2): " ++ show lat2 ++ " -> " ++ show mlat2 + + -- Get Z500 variable: get a slice, do COARDS short -> double + -- scaling, convert units. + let (Just zvar) = ncVar nc "z500" + putStrLn $ "z500 dims: " ++ show (map ncDimName $ ncVarDims zvar) + Right slice1tmp <- getA nc zvar [0, 0, 0] [1, nlat, nlon] :: SVRet CShort + let slice1 = SV.map (/ 9.8) $ coardsScale zvar slice1tmp :: SV.Vector CDouble + putStrLn $ "length slice1 = " ++ show (length slice1) + let idx i j = (nlat - i) * nlon + (j - 1) + putStrLn $ "lon(i=25) = " ++ show (lon SV.! (25 - 1)) + putStrLn $ "lat(j=40) = " ++ show (lat SV.! (nlat - 40)) + putStrLn $ "slice1(i=25,j=40) = " ++ show (slice1 SV.! idx 25 40) + +mean :: (Fractional a, Vector v a) => v a -> a +mean xs = sum xs / fromIntegral (length xs) diff --git a/examples/example2.hs b/examples/example2.hs new file mode 100644 --- /dev/null +++ b/examples/example2.hs @@ -0,0 +1,68 @@ +module Main where + +import Prelude hiding (length, maximum, minimum, sum) +import Control.Applicative ((<$>)) +import qualified Data.Map as M +import Foreign.C +import Data.Vector.Generic hiding ((++), map) +import qualified Data.Array.Repa as Repa +import Data.Array.Repa.Repr.ForeignPtr (F) +import Data.Array.Repa.Repr.Unboxed (U) +import qualified Data.Vector.Storable as SV + +import Data.NetCDF +import qualified Data.NetCDF.Vector as V +import qualified Data.NetCDF.Repa as R + +type SVRet a = IO (Either NcError (SV.Vector a)) +type FArray2 a = Repa.Array F Repa.DIM2 a +type Array2 a = Repa.Array U Repa.DIM2 a +type RepaRet2 a = IO (Either NcError (FArray2 a)) +type RepaRet1 a = IO (Either NcError (Repa.Array F Repa.DIM1 a)) + +main :: IO () +main = do + Right nc <- openFile "z500-short.nc" + + let Just ntime = ncDimLength <$> ncDim nc "time" + Just nlat = ncDimLength <$> ncDim nc "latitude" + Just nlon = ncDimLength <$> ncDim nc "longitude" + + let (Just timevar) = ncVar nc "time" + Right time <- get nc timevar :: SVRet CInt + let (Just lonvar) = ncVar nc "longitude" + Right lon <- get nc lonvar :: SVRet CFloat + let (Just latvar) = ncVar nc "latitude" + Right lat <- get nc latvar :: SVRet CFloat + let (Just zvar) = ncVar nc "z500" + + let late = vectorIndex LT FromEnd lat 40.0 + lats = vectorIndex GT FromStart lat 60.0 + lons = vectorIndex LT FromStart lon 10.0 + lone = vectorIndex GT FromEnd lon 30.0 + start = [0, lats, lons] + count = [1, late - lats + 1, lone - lons + 1] + Right slice1tmp <- getA nc zvar start count :: RepaRet2 CShort + let slice1tmp2 = coardsScale zvar slice1tmp :: FArray2 CDouble + let slice1 = (Repa.computeS $ + Repa.map (realToFrac . (/ 9.8)) slice1tmp2) :: Array2 Double + putStrLn $ "size slice1 = " ++ show (Repa.extent slice1) + putStrLn $ "maximum slice1 = " ++ show (Repa.foldAllS max (-1.0E9) slice1) + putStrLn $ "minimum slice1 = " ++ show (Repa.foldAllS min (1.0E9) slice1) + +mean :: (Fractional a, Vector v a) => v a -> a +mean xs = sum xs / fromIntegral (length xs) + +data IndexStart = FromStart | FromEnd + +vectorIndex :: (SV.Storable a, Ord a) + => Ordering -> IndexStart -> SV.Vector a -> a -> Int +vectorIndex o s v val = case (go o, s) of + (Nothing, _) -> (-1) + (Just i, FromStart) -> i + (Just i, FromEnd) -> SV.length v - 1 - i + where go LT = SV.findIndex (>= val) vord + go GT = SV.findIndex (<= val) vord + vord = case s of + FromStart -> v + FromEnd -> SV.reverse v diff --git a/examples/example3.hs b/examples/example3.hs new file mode 100644 --- /dev/null +++ b/examples/example3.hs @@ -0,0 +1,68 @@ +module Main where + +import Prelude hiding (length, sum) +import Control.Applicative ((<$>)) +import Data.NetCDF +import Data.NetCDF.HMatrix +import qualified Data.Map as M +import Foreign.C +import Foreign.Storable +import Numeric.Container + +type VRet a = IO (Either NcError (HVector a)) +type MRet a = IO (Either NcError (HRowMajorMatrix a)) + +main :: IO () +main = do + -- Open file and access some basic metadata. + Right nc <- openFile "z500-short.nc" + putStrLn $ "Name: " ++ ncName nc + putStrLn $ "Dims: " ++ show (M.keys $ ncDims nc) + putStr $ unlines $ map (\(n, s) -> " " ++ n ++ ": " ++ s) $ + M.toList $ flip M.map (ncDims nc) $ + \d -> show (ncDimLength d) ++ if ncDimUnlimited d then " (UNLIM)" else "" + putStrLn $ "Vars: " ++ show (M.keys $ ncVars nc) + putStrLn $ "Global attributes: " ++ show (M.keys $ ncAttrs nc) + + -- Get dimension sizes. + let Just ntime = ncDimLength <$> ncDim nc "time" + Just nlat = ncDimLength <$> ncDim nc "latitude" + Just nlon = ncDimLength <$> ncDim nc "longitude" + + -- Get longitude variable: metadata, then all values, then a slice + -- of values. + let (Just lonvar) = ncVar nc "longitude" + Right (HVector lon) <- get nc lonvar :: VRet CFloat + let mlon = mean lon + putStrLn $ "longitude: " ++ show lon ++ " -> " ++ show mlon + Right (HVector lon2) <- getS nc lonvar [0] [72] [2] :: VRet CFloat + let mlon2 = mean lon2 + putStrLn $ "longitude (every 2): " ++ show lon2 ++ " -> " ++ show mlon2 + + -- Get latitude variable: metadata, then all values, then a slice of + -- values. + let (Just latvar) = ncVar nc "latitude" + Right (HVector lat) <- get nc latvar :: VRet CFloat + let mlat = mean lat + putStrLn $ "latitude: " ++ show lat ++ " -> " ++ show mlat + Right (HVector lat2) <- getS nc latvar [0] [37] [2] :: VRet CFloat + let mlat2 = mean lat2 + putStrLn $ "latitude (every 2): " ++ show lat2 ++ " -> " ++ show mlat2 + + -- Get Z500 variable: get a slice, do COARDS short -> double + -- scaling, convert units. + let (Just zvar) = ncVar nc "z500" + putStrLn $ "z500 dims: " ++ show (map ncDimName $ ncVarDims zvar) + Right slice1tmp <- getA nc zvar [0, 0, 0] [1, nlat, nlon] :: MRet CShort + let (HRowMajorMatrix slice1tmp2) = + coardsScale zvar slice1tmp :: HRowMajorMatrix CDouble + slice1 = cmap ((/ 9.8) . realToFrac) slice1tmp2 :: Matrix Double + putStrLn $ "size slice1 = " ++ + show (rows slice1) ++ " x " ++ show (cols slice1) + putStrLn $ "lon(i=25) = " ++ show (lon @> (25 - 1)) + putStrLn $ "lat(j=40) = " ++ show (lat @> (nlat - 40)) + let v @!! (i, j) = v @@> (nlat - i, j - 1) + putStrLn $ "slice1(i=25,j=40) = " ++ show (slice1 @!! (25, 40)) + +mean :: (Storable a, Fractional a) => Vector a -> a +mean xs = (foldVector (+) 0 xs) / fromIntegral (dim xs) diff --git a/hnetcdf.cabal b/hnetcdf.cabal --- a/hnetcdf.cabal +++ b/hnetcdf.cabal @@ -1,5 +1,5 @@ name: hnetcdf -version: 0.2.1.0 +version: 0.2.2.0 synopsis: Haskell NetCDF library description: Bindings to the Unidata NetCDF library, along with a higher-level @@ -15,11 +15,13 @@ extra-source-files: README.markdown CHANGELOG.markdown cabal-version: >=1.8 homepage: https://github.com/ian-ross/hnetcdf +data-files: test/tst-raw-get-put.nc test/tst-bathymetry.nc source-repository head type: git location: https://github.com/ian-ross/hnetcdf + Library exposed-modules: Data.NetCDF Data.NetCDF.Vector @@ -42,7 +44,7 @@ either >= 3.4.1 && < 5, errors >= 1.4.2 && < 1.5, filepath >= 1.3 && < 1.4, - hmatrix >= 0.15.2.0 && < 0.17, + hmatrix >= 0.16.0.2 && < 0.17, repa, transformers >= 0.2 && < 0.5, vector >= 0.9 && < 0.11 @@ -108,7 +110,7 @@ base >= 4.3 && < 5, containers >= 0.4 && < 0.6, errors >= 1.4.2 && < 1.5, - hmatrix >= 0.15.2.0 && < 0.17, + hmatrix >= 0.16.0.2 && < 0.17, repa, vector >= 0.9 && < 0.11, directory >= 1.1 && < 1.3, @@ -130,3 +132,29 @@ vector >= 0.9 && < 0.11, directory >= 1.1 && < 1.3 extra-libraries: netcdf + +Executable example1 + hs-source-dirs: examples + main-is: example1.hs + build-depends: hnetcdf, + base >= 4.3 && < 5, + containers >= 0.4 && < 0.6, + vector >= 0.9 && < 0.11 + +Executable example2 + hs-source-dirs: examples + main-is: example2.hs + build-depends: hnetcdf, + base >= 4.3 && < 5, + containers >= 0.4 && < 0.6, + vector >= 0.9 && < 0.11, + repa + +Executable example3 + hs-source-dirs: examples + main-is: example3.hs + build-depends: hnetcdf, + base >= 4.3 && < 5, + containers >= 0.4 && < 0.6, + vector >= 0.9 && < 0.11, + hmatrix diff --git a/test/test-put.hs b/test/test-put.hs --- a/test/test-put.hs +++ b/test/test-put.hs @@ -14,10 +14,12 @@ main = do let xdim = NcDim "x" 10 False ydim = NcDim "y" 10 False + xattrs = M.fromList [("attr", NcAttrFloat [0.1, 0.2]), + ("name", NcAttrChar "xvar")] nc = emptyNcInfo "tst-put.nc" # addNcDim xdim # addNcDim ydim # - addNcVar (NcVar "x" NcDouble [xdim] M.empty) # + addNcVar (NcVar "x" NcDouble [xdim] xattrs) # addNcVar (NcVar "y" NcDouble [ydim] M.empty) # addNcVar (NcVar "z" NcDouble [xdim, ydim] M.empty) let write nc = do diff --git a/test/tst-bathymetry.nc b/test/tst-bathymetry.nc new file mode 100644 Binary files /dev/null and b/test/tst-bathymetry.nc differ diff --git a/test/tst-raw-get-put.nc b/test/tst-raw-get-put.nc new file mode 100644 Binary files /dev/null and b/test/tst-raw-get-put.nc differ