MATLAB Answers

0

NETCDF GIS coupling creationg problems

Asked by Coastal_Eng on 29 Oct 2013
Latest activity Commented on by Coastal_Eng on 4 Nov 2013
Hi guys. I'm having problems in creating a proper netcdf file in matlab. I'm able to create the file and the time-series, but i have problems with the axis, and when I try to project my netCDF it has a very weird aspect.
So let me tell you about the input. I have the X and Y coord already in grid format:
( X(72x888) Y(72x888)),
I have the time
endT=(49x1)
and finally i have my time series:
value=49x72x888
my objective is to build a netcdf in such a way that i can easily export the time serie to arcGIS and build an animation. The problem i'm facing is that when exported to gis the projection is messed up, here is how it looks:
and here is how it should look:
I think my code is messed up on the axis genereation, but i'm unable to detect the problem, could you guys give me a hand?! Here is the code:
clc; clear all; close all;
load('water level.mat')
%organizing and visualizing
X=data.X;
Y=data.Y;
T=data.Time;
value=data.Val;
X=X(2:end,2:end);
Y=Y(2:end,2:end);
w1=value(1,:,:);
w1=squeeze(w1);
xyz=[X(:) Y(:) w1(:)];
plot(xyz(:,1),xyz(:,2),'.')
axis equal
pcolor(X,Y,w1)
shading interp
axis equal
%starting with the netCDF creationg
scope = netcdf.create('WATER_LEVEL_final.nc','netcdf4');
NC_GLOBAL = netcdf.getConstant('NC_GLOBAL');
fillValue = -9999;
dimidX1 = netcdf.defDim(scope,'x1_coor',length(w1(:,1)));
dimidY1 = netcdf.defDim(scope,'y2_coor',length(w1));
dimidT = netcdf.defDim(scope,'time',length(T));
netcdf.putAtt(scope,NC_GLOBAL,'title','FILL LATER')
netcdf.putAtt(scope,NC_GLOBAL,'long_title','FILL LATER')
netcdf.putAtt(scope,NC_GLOBAL,'comments','FILL LATER')
netcdf.putAtt(scope,NC_GLOBAL,'institution','LPO/Ifremer')
netcdf.putAtt(scope,NC_GLOBAL,'source','FILL LATER')
netcdf.putAtt(scope,NC_GLOBAL,'LATER','LATER')
netcdf.putAtt(scope,NC_GLOBAL,'LATER','LATER')
netcdf.putAtt(scope,NC_GLOBAL,'CreationDate',datestr(now,'yyyy/mm/dd HH:MM:SS'))
netcdf.putAtt(scope,NC_GLOBAL,'CreatedBy',getenv('RODRIGO GONÇALVES'))
%Then add the variable:
varid = netcdf.defVar(scope,'time','double',[dimidT]);
netcdf.putAtt(scope,varid,'standard_name','time');
netcdf.putAtt(scope,varid,'long_name','Time of measurements');
netcdf.putAtt(scope,varid,'units','FILL LATER ');
netcdf.defVarFill(scope,varid,false,fillValue);
netcdf.putVar(scope,varid,T);
D = value;
varname = 'value_water_level';
long_name = 'This is a time serie of noise !';
unit = 'm';
% I THINK PROBLEM IS HERE !
varid = netcdf.defVar(scope,varname,'double',[dimidT,dimidX1,dimidY1]);
netcdf.putAtt(scope,varid,'long_name',long_name);
netcdf.putAtt(scope,varid,'units',unit);
netcdf.defVarFill(scope,varid,false,fillValue);
v = D;
v(isnan(v)) = fillValue;
netcdf.putVar(scope,varid,v);
netcdf.close(scope)
Tanks ! ROD

  0 Comments

Sign in to comment.

1 Answer

Answer by Ashish Uthama on 29 Oct 2013

I am not familiar with arcGIS, does it show you the variable dimensions? Does it expect time to be first or the last dimension?
As you guess, its either the order of the dimesion IDs in the defVar call, or the order of the data in v. If arcGIS expects the dimension ordering of [X, Y, T], you should use [Y, X, T] in MATLAB. (MATLAB is colum major, i.e the first dimension is the column (Y), other programs based on C are usually row major (X is the first dimension). You could use reshape to change the order of data in v to match your dimension ordering.
Maybe you could try with smaller data like v = reshape(1:5*4*3,[5 4 3]) to better puzzle this out?

  3 Comments

Hey, Thaks for your answer. arcGIS is able to directly load the netcdf provided from NOAA for axample. are you familiar with those files?
I can also tell you, that if i plot a squeezed dimension of the original 49x72x888 variable i get exactly what is shown in the GIS.
X=data.X;
Y=data.Y;
T=data.Time;
value=data.Val;
X=X(2:end,2:end);
Y=Y(2:end,2:end);
w1=value(1,:,:);
w1=squeeze(w1);
xyz=[X(:) Y(:) w1(:)];
plot(xyz(:,1),xyz(:,2),'.')
axis equal
pcolor(w1)
shading interp
axis equal
This makes me think that somehow i must insert the original coordinate matrices into the netcdf. In sintax my time serie variable, should be X Y VALUE, where X=(72x888) Y=(72x888) and value=49x72x888. But i'm not able to do that. Do you understand my doubt !? cheers
I guess you need to create a variable with the same name as the dimension to store the coordinates, i.e a Coordinate variable
Did that, but seems that then i got a problem with the x and y coord, spacing. they must be the same in order to import it in gis. I dnt know what to do anymore. Tks for your assistance.

Sign in to comment.