[Scilab-users] Getting netCDF files into Scilab 6.0.1

[arctica1963]

Hi all,

Just a quick update. Used a function to import an ESRI ascii raster seems to
work ok:

function [X,Y,Z]=asc_2_xyz(filename)
    
    fd=mopen(filename,'r'); // open ascii raster for read only
    k=mfscanf(fd,'%s\t%d'); nc=k(2); // number of columns
    k=mfscanf(fd,'%s\t%d'); nr=k(2); // number of rows
    k=mfscanf(fd,'%s\t%f'); x1=k(2); // X coordinate = lower left
    k=mfscanf(fd,'%s\t%f'); y1=k(2); // Y coordinate = lower left
    k=mfscanf(fd,'%s\t%f'); dxy=k(2); // X and Y increment
    k=mfscanf(fd,'%s\t%f'); nodata=k(2); // No data value (e.g. -9999)
    mclose(fd)

// x1 and y1 = minimum X and Y

    x2 = x1+dxy*nc; // Define maximum X coordinate
    y2 = y1+dxy*nr; // Define maximum Y coordinate

    Z=fscanfMat(filename); // Z matrix
    Z(Z==nodata)=%nan; // if any 'nodata' values in grid set to internal NaN
value
    Z=flipdim(Z,1); // for some reason the data are flipped on Y on import
so needed flipping

    xx=linspace(x1,x2,nc); // same number of columns but increment different
from dxy
    yy=linspace(y1,y2,nr);
    [X,Y]=meshgrid(xx,yy); // Coordinate matrices
    
endfunction

All looks fine, but the increment from linspace cannot match dxy unless one
has nc+1. Not a big issue as this can be corrected in GMT. But also reversed
the thing to output .asc and that gets all the dimensions and increment
fine:

function xyz_2_asc(filename,X,Y,Z,nodata,n)
// n = precision (e.g. 6)
// nodata should be appropriate for the Z data range
//e.g. xyz_2_asc(test.asc,x,y,z,-9999,6)
    
// Define matrix dimensions
    dimx=size(X); // 2-D
    dimy=size(Y);
    dimz=size(Z);
    
// Grid properties
    minx=min(X); //initial X coordinate (Lower left)
    maxx=max(X);
    miny=min(Y); //initial Y coordinate (Lower left)
    maxy=max(Y);
    
    dxy=abs(maxx-minx)/dimx(1) // correct for increment precision issues
    
// Z matrix
    Z(isnan(Z))=nodata; // add nan values (where present)) to output file
    Z=flipdim(Z,1); // invert matrix rows in Y
    
// Output file
    
    fd=mopen(filename,'w'); // open file for write
// Create ESRI Ascii raster header
//
// ncols 120
// nrows 120
// xllcorner -5
// yllcorner -5
// cellsize 0.0166666666667 - must be square pixels for ESRI ascii raster
// nodata_value -9999 (select value appropriate for dataset)
//
    mfprintf(fd,'ncols\t');
    mfprintf(fd,'%d\n',dimz(2));
    mfprintf(fd,'nrows\t');
    mfprintf(fd,'%d\n',dimz(1));
    mfprintf(fd,'xllcorner\t');
    mfprintf(fd,'%f\n',minx);
    mfprintf(fd,'yllcorner\t');
    mfprintf(fd,'%f\n',miny);
    mfprintf(fd,'cellsize\t');
    mfprintf(fd,'%.10f\n',dxy); // rounding errors > %.10f but not critical
to output
    mfprintf(fd,'NODATA_value\t');
    mfprintf(fd,%f\n',nodata);

// write Z values to file	
	
    for i=1:dimz(1)
        for j=1:dimz(2)
            if j==dimz(2)
                mfprintf(fd,'%f\n',Z(i,j));
            else
                mfprintf(fd,'%f\t',Z(i,j));
            end
        end
    end
    mclose(fd);
	
endfunction

Ideally, it would be good to work out a way of reading the binary netCDF
grids, much smaller size, but at least it seems to work albeit a bit quirky
with the increments and flipping!

Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html