将大型 xyz 文件转换为网格化数据 (Matlab)
Convert large xyz file into gridded data (Matlab)
我有一个很大的 XYZ 文件(300276x3,这个文件包括 x 和 y 坐标(不是 lat/lon,而是极地立体)和高程 z),我想知道是否可以转换它进入网格数据集(n x m 矩阵)。 xyz 文件可以从以下网址下载:
并通过以下方式导入 matlab:
AIS_SEC = importdata('AIS_SEC.xyz');
我试过:
X= XYZ(:,1);
Y= XYZ(:,2);
Z= XYZ(:,3);
xr = sort(unique(X));
yr = sort(unique(Y));
gRho = zeros(length(yr),length(xr));
gRho = griddata(X,Y,Z,xr,yr')
imagesc(gRho)
Requested 300276x300276 (671.8GB) array exceeds maximum array size preference. Creation of arrays
greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size
limit or preference panel for more information.
我试过了:
% Get coordinate vectors
x = unique(XYZ(:,1)) ;
y = unique(XYZ(:,2)) ;
% dimensions of the data
nx = length(x) ;
ny = length(y) ;
% Frame matrix of grid
D = reshape(XYZ(:,3),[ny,nx]) ;
% flip matrix to adjust for plot
H = flipud(H) ;
% Transpose the matrix
H = H' ; % Check if is required
surf(x,y,H) ;
Error using reshape
To RESHAPE the number of elements must not change.
我现在可以使用 scatter3 绘制 nx3 文件(见图)
scatter3(XYZ(:,1),XYZ(:,2),XYZ(:,3),2,XYZ(:,3)) ;
colorbar
但我想用 imagesc 来做。因此,我想将 nx3 文件转换为 nxm 矩阵(raster/gridded 格式),另外我希望它作为 geotiff 文件用于 QGIS。
谢谢!
您快到了...查看您收到的有关数组大小的消息,unique(X)
的结果似乎可能会产生 300276 个唯一值,这可能是由于一些嘈杂的数据。
因此,您可以在需要的域上定义一些新向量,而不是将 griddata
与这些大型 X
和 Y
向量一起使用:
% make some sample data
N = 1000;
xv = linspace(-10,10,N);
yv = linspace(-10,10,N);
[XV,YV] = meshgrid(xv,yv);
ZV = XV.^2 + YV.^2;
% make into long vectors:
X = XV(:);
Y = YV(:);
Z = ZV(:);
% make x and y vector to interpolate z
N = 50; % size of new grid
xv = linspace(min(X), max(X), N);
yv = linspace(min(Y), max(Y), N);
[XV,YV] = meshgrid(xv,yv);
% use griddata to find right Z for each x,y pair
ZV_grid = griddata(X,Y,Z,XV,YV);
% look at result
figure();
subplot(211)
imagesc(ZV);
subplot(212);
imagesc(ZV_grid)
我有一个很大的 XYZ 文件(300276x3,这个文件包括 x 和 y 坐标(不是 lat/lon,而是极地立体)和高程 z),我想知道是否可以转换它进入网格数据集(n x m 矩阵)。 xyz 文件可以从以下网址下载:
并通过以下方式导入 matlab:
AIS_SEC = importdata('AIS_SEC.xyz');
我试过:
X= XYZ(:,1);
Y= XYZ(:,2);
Z= XYZ(:,3);
xr = sort(unique(X));
yr = sort(unique(Y));
gRho = zeros(length(yr),length(xr));
gRho = griddata(X,Y,Z,xr,yr')
imagesc(gRho)
Requested 300276x300276 (671.8GB) array exceeds maximum array size preference. Creation of arrays
greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size
limit or preference panel for more information.
我试过了:
% Get coordinate vectors
x = unique(XYZ(:,1)) ;
y = unique(XYZ(:,2)) ;
% dimensions of the data
nx = length(x) ;
ny = length(y) ;
% Frame matrix of grid
D = reshape(XYZ(:,3),[ny,nx]) ;
% flip matrix to adjust for plot
H = flipud(H) ;
% Transpose the matrix
H = H' ; % Check if is required
surf(x,y,H) ;
Error using reshape
To RESHAPE the number of elements must not change.
我现在可以使用 scatter3 绘制 nx3 文件(见图)
scatter3(XYZ(:,1),XYZ(:,2),XYZ(:,3),2,XYZ(:,3)) ;
colorbar
但我想用 imagesc 来做。因此,我想将 nx3 文件转换为 nxm 矩阵(raster/gridded 格式),另外我希望它作为 geotiff 文件用于 QGIS。
谢谢!
您快到了...查看您收到的有关数组大小的消息,unique(X)
的结果似乎可能会产生 300276 个唯一值,这可能是由于一些嘈杂的数据。
因此,您可以在需要的域上定义一些新向量,而不是将 griddata
与这些大型 X
和 Y
向量一起使用:
% make some sample data
N = 1000;
xv = linspace(-10,10,N);
yv = linspace(-10,10,N);
[XV,YV] = meshgrid(xv,yv);
ZV = XV.^2 + YV.^2;
% make into long vectors:
X = XV(:);
Y = YV(:);
Z = ZV(:);
% make x and y vector to interpolate z
N = 50; % size of new grid
xv = linspace(min(X), max(X), N);
yv = linspace(min(Y), max(Y), N);
[XV,YV] = meshgrid(xv,yv);
% use griddata to find right Z for each x,y pair
ZV_grid = griddata(X,Y,Z,XV,YV);
% look at result
figure();
subplot(211)
imagesc(ZV);
subplot(212);
imagesc(ZV_grid)