PG2PolandGeoid.m

clc

egmdatnum = importdata('EGM96PolandGeoid.gdf',' ',36);
egmdat = egmdatnum.textdata;
egmnum = egmdatnum.data;
la=egmnum(:,1);
fi=egmnum(:,2);
dat1=egmnum(:,3);
lamin=min(la);
lamax=max(la);
fimin=min(fi);
fimax=max(fi);

##plot(la,fi,dat1,'Color',[.4 .4 .4],'LineWidth',1);
##hold on

la1=la(1);la2=la(2);
  i=0;
  while la1==la2 
    i=i+1; 
    la2=la(2+i);
  end 
  fi1=fi(1);fi2=fi(2);
  j=0;
  while fi1==fi2
    j=j+1;
    fi2=fi(2+j);
  end
  dla=abs(la1-la2);
  dfi=abs(fi1-fi2);
  frows=abs(floor((fimax-fimin)/dfi+1.5));
  lcols=abs(floor((lamax-lamin)/dla+1.5));
  n=1;
  for i=frows:-1:1
    for j=1:lcols
      if n <= length(la) 
        if la(n) < 180
         X(i,j)=la(n);
        else
          X(i,j)=la(n)-360;
        end
      Y(i,j)=fi(n);
      Z(i,j) = dat1(n);
      n=n+1;
      end
    end    
  end
  
[c,h]=contourf(X,Y,Z,1.9409290E+01:1.22427358:5.0761796E+01);
clabel(c,h);
hold on

colorbar

ylabel('latitude')
xlabel('longitude')

bfid = fopen('polandn.bln');

while feof(bfid)==0;
    
    n = fscanf(bfid,'%i,%i',2);
    if isempty(n)==1;
        break
    end
    fila = fscanf(bfid,'%g,%g',[2 n(1)]);
    if n(1)~=length(fila(1,:));
        fila = fscanf(bfid,'%g %g',[2 n(1)]);
    end
    la=fila(1,:);
    fi=fila(2,:);
    plot(la,fi,'Color',[.4 .4 .4],'LineWidth',1);
    hold on

end