Skip to content

Commit

Permalink
reading of GeoSoft files
Browse files Browse the repository at this point in the history
  • Loading branch information
@cultpenguin
cultpenguin committed Apr 17, 2020
1 parent 09c0125 commit 7e18ba1
Show file tree
Hide file tree
Showing 3 changed files with 279 additions and 0 deletions.
57 changes: 57 additions & 0 deletions geosoft_xyz_plot.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
% geosoft_xyz_plot(D,LINE,BOT,TOP,ilines_plot);
%
% Call:
% geosoft_xyz_plot(D,LINE,BOT,TOP,ilines_plot);
%
% See also: geosoft_xyz_read
%
function geosoft_xyz_plot(D,LINE,BOT,TOP,ilines_plot);
if nargin<5
ilines_plot=1;
end

clim=[1 3];
cmap=colormap;
ncolors=size(cmap,1);


nlayers=min([length(BOT) size(D{1}.RES,1)]);

ids=20;
lw=3;

for iline=ilines_plot;
figure(iline);clf;set_paper;

nlayers=min([length(BOT) size(D{1}.RES,1)]);

ns=size(D{iline}.RES,2);
%imagesc(X,BOT,log10(D{iline}.RES));

% loop over
for is=1:ids:ns;
X=is;
%X=D{iline}.UTMX(is);
%X=D{iline}.UTMY(is);
for il=1:nlayers;
try
logres=log10(D{iline}.RES(il,is));
if ~isnan(logres)
icol=interp1(linspace(clim(1),clim(2),ncolors),1:ncolors,logres,'nearest','extrap');
col=cmap(icol,:);
plot([1 1].*X,D{iline}.ELEV(is)+[BOT(il) TOP(il)],'k-','color',col,'LineWidth',lw);
if (is==1&il==1);hold on;end
end
catch
end
end

end
hold off
daspect([10 1 1])
txt=sprintf('Line%d',LINE(iline));
title(txt)
try
print_mul(txt)
end
end
151 changes: 151 additions & 0 deletions geosoft_xyz_read.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
% geosoft_xyz_read: reads Geosoft XYZ ASCII data format
% Call
%
% % column header in line 2, and no strings in data section
% [D,LINE,CHEAD,HEADER]=geosoft_xyz_read(filename);
%
% % column header in line 5, and no strings in data section
% chead_line=5;
% [D,LINE,CHEAD,HEADER]=geosoft_xyz_read(filename,chead_line);
%
% % column header in line 5, and possible strings in data section
% chead_line=5;
% lookfor_strings=1; %(SLOW)
% [D,LINE,CHEAD,HEADER]=geosoft_xyz_read(filename,chead_line,no_strings);
%
%
%
function [D,LINE,CHEAD,HEADER]=geosoft_xyz_read(filename,chead_line,lookfor_strings);
if nargin<1;
filename='YukonAEM_INV.XYZ';
end
if nargin<2;
chead_line=2;
disp(sprintf('%s: Assuming column header in line %d',mfilename,chead_line));
end
if nargin<3;
lookfor_strings=0;
disp(sprintf('%s: Assuming there are no strings in data!!',mfilename));
end


CHEAD=[];
disp(sprintf('%s: Reading %s',mfilename,filename))
fid=fopen(filename,'r');


NC=100;
NL=28;
IC_X=1;
IC_Y=2;
IC_ELEV=3;
IC_SKIP=4;

IC_BOT=IC_SKIP+[1:NL];
IC_TOP=IC_SKIP+NL+[1:NL];
IC_RES=IC_SKIP+2*NL+[1:NL];

%USE_LINES=10142;
USE_LINES='';;
USELINE=0;
iline=0;
iheader=0;
line_in_file=0;
while ~feof(fid);
line_in_file=line_in_file+1;
l=fgetl(fid);
if strcmp(lower(l(1)),'/');
% COMMENT OR HEADER
is_header=0;
if length(l)==1;
is_header=1;
else
if strcmp(lower(l(2)),'/')
is_header=0; % is_comment=1;
else
is_header=1;
end
end
if is_header==1
disp(sprintf('%s: Reading header line %d: ',mfilename,iheader+1,l));
iheader=iheader+1;
HEADER{iheader}=l;
if (iheader==chead_line);
CHEAD=strsplit(strtrim(l(2:end)));
end
else
disp(sprintf('%s: Skipping line %d: ',mfilename,i,l));
end

elseif (~strcmp(l(1),' '));%(isempty(str2num(l(1)))); %strcmp(lower(l(1:4)),'line');

LINENUMBER=str2num(l(5:end));
if line_in_file==12, keyboard;end
if ((~isempty(find(LINENUMBER==USE_LINES)))|(length(USE_LINES)==0));
USELINE=1;
disp(sprintf('reading line : %d',LINENUMBER));;
iline=iline+1;
LINE(iline)=LINENUMBER;

else
USELINE=0;
disp(sprintf('skipping line : %d',LINENUMBER));;
end

idata=0;

elseif (USELINE==1);
% ACTUAL DATA
idata=idata+1;
if lookfor_strings==0;
% assume no strings in data section
linedata=sscanf(l,'%f');
else
% possible strings in data sections (much slower)
line_split=strsplit(strtrim(l));
for k=1:length(line_split);
linedata(k)=str2num(line_split{k});
end
end

% if idata==1;
% BOT=linedata(IC_BOT);
% TOP=linedata(IC_TOP);
% end
% n_linedata=length(linedata);
% RES=linedata(IC_RES(1):end);
% D{iline}.RES(1:NL,idata)=NaN;;
% D{iline}.RES(1:length(RES),idata)=RES;
% D{iline}.ELEV(idata)=linedata(IC_ELEV);
% D{iline}.UTMX(idata)=linedata(IC_X);
% D{iline}.UTMY(idata)=linedata(IC_Y);
% %D{iline}(idata,1:n_linedata)=sscanf(l,'%f');

try
D{iline}(idata,1:length(linedata))=linedata;
catch
keyboard
end
else
keyboard
% do nothin
end


end


fclose(fid);
try
if nargout==0;
[p,f,e]=fileparts(filename);
if (length(USE_LINES)==1)
f_out=sprintf('%s_%g',f,USE_LINES);
else
f_out=[f,'.mat'];
end
save(f_out,'D','LINE','TOP','BOT','H');
end
catch
keyboard;
end
71 changes: 71 additions & 0 deletions geosoft_xyz_scatter.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
% geosoft_xyz_scatter
%
% Call:
% geosoft_xyz_scatter(D,id_array,ix,iy,CHEAD,ilines_plot,do_log);
%
% See also: geosoft_xyz_read
%
function geosoft_xyz_scatter(D,id_array,ix,iy,CHEAD,ilines_plot,do_log);
if nargin<3, ix=1; end
if nargin<4, iy=2; end
if nargin<2, id_array=[]; end
if nargin<6
ilines_plot=1:length(D);
end
if isempty(ilines_plot);ilines_plot=1:length(D);end
if nargin<6, do_log=1;end



if isempty(id_array);
for i=ilines_plot;
i;
plot(D{i}(:,ix),D{i}(:,iy),'-');
hold on
end
hold off
return
end


ssize=2;


nsp=ceil(sqrt(length(id_array)));

j=0;
for id=id_array
j=j+1;

if length(id_array)>1
subplot(nsp,nsp,j);
end



for iline=ilines_plot

i_ok=find( (D{iline}(:,ix)~=0) & (D{iline}(:,ix)~=0) );
i_max=3000;
ii=i_ok(unique(ceil(linspace(1,length(i_ok),i_max))));



if (do_log==1);
scatter(D{iline}(ii,ix), D{iline}(ii,iy),ssize,log(real(D{iline}(ii,id))),'filled')
else
scatter(D{iline}(ii,ix), D{iline}(ii,iy),ssize,(real(D{iline}(ii,id))),'filled')
end
hold on
end
hold off
axis image
xlabel(CHEAD{ix})
ylabel(CHEAD{iy})
if exist('CHEAD','var')
title(sprintf('%s (icol=%d)',CHEAD{id},id),'Interp','none')
end
colormap(cmap_geosoft)
colorbar
drawnow;
end

0 comments on commit 7e18ba1

Please sign in to comment.