-
Notifications
You must be signed in to change notification settings - Fork 1
/
PG2TurkeyGeoid.m
75 lines (65 loc) · 1.3 KB
/
PG2TurkeyGeoid.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
clc
egmdatnum = importdata('EGM96TurkeyGeoid.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.0607855E+02:1.90811467:8.4732917E+01);
clabel(c,h);
hold on
colorbar
ylabel('latitude')
xlabel('longitude')
bfid = fopen('turkey.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