Convert Matrix to ArcGIS raster

Hello Everybody,

I am working in a project, where a team is generating simulations with a model and exporting results in .mat files (x, y, z). I need to plot the results with ArcGIS (only ploting without doing further interpolations or operations, just setting colors and legend).

The problem is I need to generate a raster in matlab so ArcGIS can read it. If I only import x,y,z points into ArcGIS I will need to create a shapefile with them and interpolate it to generate a raster, which is not the best option, as it can differ from Matlab results (interpolation method etc).

So, anyone could help me converting a mx3 matlab variable (containing x, y coordinates in WGS84) and a variable z into a raster so I can read it directly into ArcGIS as a raster file?

Thanks in advance.

Best regards,

Diego

Answer 1

[KSSV]

It depends whether your data is structured or unstructured. Based on the data, the folloiwng is the procedure. Let A be your (x,y,z) data.

x = A(:,1) ; y = A(:,2) ; z = A(:,3) ;
%%structured 
xi = unique(x) ; yi = unique(y) ;
[X,Y] = meshgrid(xi,yi) ;
Z = reshape(z,size(X)) ;
figure
surf(X,Y,Z)
%%unstructured 
Z = griddata(x,y,z,X,Y) ; 

[Diego Fonseca]

Thanks for your answer, KSSV.

However, I have aleady a griddded data. The question is how to convert this dataset (my mx3 x,y,z, where x and y are WGS84 coordinates) into a file, I can read directly as a raster into ArcMap (I mean, I will not need to interpolate the data with ArcGIS to create my own raster).

[KSSV]

ArcMap reads netCDF file right, you can write the raster to netCDF and use it in ArcMap.

[Diego Fonseca]

I have never done that. Do you know how i can perform this procedure in MAtlab from my x,y,z data matrix to a netCDF file?

[KSSV]

% nc filename to be written
file = 'myfile.nc' ;

%% Write lon and lat variables
% Get data
lon = xi ;
lat = yi ;
nx = length(lon) ;
nccreate(file,'lon','Dimensions',{'lon',1,nx},'DeflateLevel',7) ;
ny = length(lat) ;
nccreate(file,'lat','Dimensions',{'lat',1,ny},'DeflateLevel',7) ;
nccreate(file,'time','Dimensions',{'time',1,Inf},'DeflateLevel',7) ;
nccreate(file,'z','Dimensions',{'lon','lat'},'DeflateLevel',7) ; 
ncwrite(file,'z',Z,[1,1]) ;   % write 2D data 

[Diego Fonseca]

Hi KSSV. Thanks fot the answer and sorry for the time to accept it.

I was able to create a netCDF file and read it in ArcGIS. There, I coud project it.

Only problem is that the code you indicated created an extra point in the netCDF file, with a z value of 9,96921e+036.

Any suggestion how to avoid it? Moreover, do you have any idea if I can directly define the coordinate system to WGS84 in matlab, in order to speed up the process (instead of importing it to ArcGIS and defineing the coordinate system there)?

[KSSV]

WGS84 means UTM coordinates?

Just use

ncwrite(file,'z',Z) ;   % write 2D data 

[Diego Fonseca]

No, UTM is the projection, WGS 84 is the datum for geographic coordinates. However, this is not a major issue, as I can project the netcdf in arcgis (it was just to speed up the process).

Answer 2

[Robert Daly]

Following on from KSSV in the comments...

You can export the data to a netCDF file then import that in ArcGIS.

To define the coordinate system for ArcGIS you need to add a bunch of attributes to the variables.

To work out the important attributes I exported a source raster to netCDF using "Multidimension tools>Raster to NetCDF"

Then in matlab use

ncdisp('export.nc')

to read the attributes.

The "esri_pe_string" attribute value is important according to this article Stackexchange

I got the string from the exported example.

% nc filename to be written
file = 'test.nc' ;
% create test dataset
Z=rand(10,10);
x = [2.2493e+05:2.2493e+05+10]
y = [6.1817e+06:6.1817e+06+10]
[x,y]=meshgrid(x,y)
% Write X and Y variables
% Get data
X = x(1,:) ;
Y = y(:,1) ;

%% Define data dimensions and write data
nccreate(file,'x','Dimensions',{'x',1,length(X)},'DeflateLevel',5) ;
ncwrite(file,'x',X) ;
nccreate(file,'y','Dimensions',{'y',1,length(Y)},'DeflateLevel',5) ;
ncwrite(file,'y',Y) ;
nccreate(file,'z','Dimensions',{'x','y'},'DeflateLevel',5) ; 
ncwrite(file,'z',Z',[1,1]) ;   % write 2D data 
nccreate(file,'transverse_mercator','Dimensions',{'transverse_mercator',1,Inf},'DeflateLevel',5) ;

% Add projection attributes for ArcGIS
ncid = netcdf.open(file,'WRITE')
netcdf.reDef(ncid);

netcdf.putAtt(ncid,0,'units','Meter')
netcdf.putAtt(ncid,0,'long_name', 'Easting')
netcdf.putAtt(ncid,0,'standard_name','projection_x_coordinate')
netcdf.putAtt(ncid,0,'axis','X')

netcdf.putAtt(ncid,1,'units','Meter')
netcdf.putAtt(ncid,1,'long_name', 'Northing')
netcdf.putAtt(ncid,1,'standard_name','projection_y_coordinate')
netcdf.putAtt(ncid,1,'axis','Y')

netcdf.putAtt(ncid,2,'long_name','Habitat Suitability Index')
netcdf.putAtt(ncid,2,'standard_name','HSI')
netcdf.putAtt(ncid,2,'esri_pe_string','PROJCS["GDA2020_MGA_Zone_54",GEOGCS["GDA2020",DATUM["GDA2020",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",141.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]')
netcdf.putAtt(ncid,2,'coordinates','x y')
netcdf.putAtt(ncid,2,'grid_mapping','transverse_mercator')

netcdf.putAtt(ncid,3,'grid_mapping_name','transverse_mercator')
netcdf.putAtt(ncid,3,'longitude_of_central_meridian',141)
netcdf.putAtt(ncid,3,'latitude_of_projection_origin',0)
netcdf.putAtt(ncid,3,'scale_factor_at_central_meridian',0.9996)
netcdf.putAtt(ncid,3,'false_easting',500000)
netcdf.putAtt(ncid,3,'false_northing',10000000)

netcdf.close(ncid)