SAS proc gproject et probleme de projection [Résolu]

[JBLamy]

Bonjour à tous,

J'ai un jeu de donnée de coordonnées géographiques (en degree decimal) mise au point sur un logiciel type SIG (GIS software (TNTmips©, Microimage LTD)) que j'appelle "ob_anno" dans la suite. Ces points sont située en europe continental.

Quand j'affiche "ob_anno" sur googleearth, pas de probleme (voir fig 1) ! Les points sont bien là ou ils devraient être. Google earth n'affiche pas tout les points car ils y en a 400 dans le fichiers et il ne peut qu'en afficher que 200 mais j'ai verifier que tout les points étaient correctement projecté sur google earth.
En revanche, quand je projette ces points sur les cartes fournis par SAS (voir fig 2) ... Les points sont décalés certains se retrouvent dans la mer... et ce quelque soit la projection utilisée.

Voici la partie du code incriminée

Code :

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
*Import the world map AND convert the radian IN degree ;
 
* je convertis les coordonnées en radians de la carte SAS en degree afin de fusionner les deux jeu de données, mes points (en degréé aussi) + la carte sas pour faire la projection sur le même jeu de donnée (ça reduit la distortion de faire la projection en une fois) ; 
 
DATA work.world ; 
SET maps.world ;
x = long/(constant('pi')/180) ;
label x = x ;  
y = lat/(constant('pi')/180) ;
label y = y ;
run ; 
 
*SET world map + pinus pinaster DATA TO get a good projection ;
 DATA ALL ; 
 length id 8 ; 
 SET world ob_anno ; 
 run ; 
 
proc sort DATA = allr ; 
BY var x y ; 
run ; 
 
%let pro = ortho ;
proc gproject DATA= ALL out = all_p degree latmin = 30  latmax = 55   longmin = -15  longmax = 10 project = &pro dupok  ; 
id id ;
run ; 
 
DATA world1p ob_mp ; 
SET all_p ;
IF id = -999 then output ob_mp ; 
else output world1p ; 
run ; 
 
proc sort DATA = world1p ; BY id; run ; 
 
proc gmap DATA = world1p map = world1p ;
id id ;
choro id / nolegend levels=1  coutline=graycc anno= ob_mp  ;
run;
quit;

[datametric]

Peux-tu

nous donner un bout de tes données pour tester
vérifier pourquoi tu fais une proc sort sur ALLR et non sur ALL
Recharger les images car aucune pièce n'est jointe.

[JBLamy]

Bonjour,

--Premierement, merci Datametric de me consacrer un peu de votre temps...

--Secondement, le "proc sort" sur "allr" peut être remplacer par "all", c'est simplement un copier coller un peu rapide, je faisais des tests sur un jeu de donnée réduit (d'ou le r). Donc :

proc sort data = allr ; => proc sort data = all ;

--Les figures sont rechargées... fig 1 = "google_earth.pdf" ; fig 2 = Pinus_pinaster_back_groundJPEG.jpeg (attention on ne voit pas bien les points noires correspondant à la projection, c'est un autre détail que je dois regler).
--L'input ok , demain, si cela ne marche pas, j'en reverrais une autre version parce que je n'ai pas SAS sur cet ordi... j'ai reconstruis l'input au plus proche de l'input correspondant au code fournis.
à bientôt Datametric !

JB

[datametric]

j'obtiens une carte toute griffonnée

tu peux monter ta table ob_anno ?

[JBLamy]

Envoyé par datametric

j'obtiens une carte toute griffonnée

tu peux monter ta table ob_anno ?

J'ai eu le même soucis ... en fait cela ne vient pas de la table anno mais plutôt des "drop" ou des "proc sort" que l'on fait sur les cartes de SAS. je t'envoie le code entier de cette partie là (voir "prog.txt").
J'ajoute la table "ob_anno" entiere (voir "esSSSSai.txt")...
A de suite, en esperant que t'y verrait plmus clair.

[JBLamy]

Je vois que personne ne réagit, c'est peut être que les informations sont trop disparates. Donc ici :
Le programme utilisé et testé sur SAS 9.2, il est destiné à mettre des points sur une carte d'europe. Les points (en coordonnées géographiques déciamle) a projeter ont été projeté avec google earth, ils apparaissent bien ou ils doivent être (cf image : google_earth_projection.pdf).
En revanche, cela ne marche pas avec le code (cf Prog_map) car cela donne l'image suivante ( Pinus_pinaster_back_groundJPEG.JPG). Je fournis la table ob_anno qui sont l'ensemble des points à projeter.
Merci,

[JBLamy]

merci a Robert Allison pour l'aide sur ce code :

Code :

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
DATA world; SET maps.world;
label x='x';
label y='y';
x=long;
y=lat;
run ;
 
DATA ob_anno;
input lat long;
/* convert lat/long from degrees to radians, so they are same as map */
/* also, multiply longitude by -1, so it is 'westlong' like the map */
x=atan(1)/45*long*-1;
y=atan(1)/45*lat;
flag=1;
datalines;
32.1404 -5.44198
32.4662 -5.05228
32.4662 -5.43449
32.4662 -5.78729
32.4913 -4.28729
32.4913 -4.69890
32.8422 -3.86762
32.8422 -4.27922
32.8422 -4.66143
32.8422 -5.04364
32.8422 -5.42584
33.1932 -4.24175
33.1932 -5.03556
33.2182 -4.62338
33.5441 -4.61589
33.5441 -5.38030
33.5441 -6.14471
33.5692 -3.88030
33.5692 -4.23311
33.5692 -4.99752
33.5943 -3.52692
33.8951 -3.84341
33.8951 -4.22561
33.8951 -4.57842
33.8951 -4.99002
33.8951 -6.16604
33.9201 -3.49002
33.9201 -5.37165
33.9201 -5.75386
34.2711 -4.62857
34.2711 -4.98138
34.2961 -3.80478
34.2961 -5.77461
34.5970 -3.47446
34.5970 -4.20947
34.5970 -4.62108
34.6220 -3.85609
34.9479 -3.46639
34.9479 -5.34802
34.9730 -3.84802
34.9730 -4.25962
34.9730 -4.55363
34.9730 -4.93583
35.2738 -5.39932
35.3239 -4.95716
35.6498 -4.94967
35.6498 -5.33188
36.3767 6.41563
36.3767 6.00402
36.3767 -4.90355
36.3767 -5.31516
36.4018 4.50517
36.4018 4.12297
36.4018 3.32915
36.4269 3.74134
36.7277 6.77650
36.7277 6.01209
36.7277 5.62988
36.7277 3.74826
36.7277 -4.57208
36.7277 -4.89548
36.7527 6.45367
36.7527 4.54264
36.7527 4.10164
36.7527 -5.33591
36.8280 8.77804
36.9032 -4.30344
37.0786 6.81397
37.0786 -4.18180
37.1037 6.46174
37.1037 -3.03460
37.1037 -3.44621
37.1037 -3.79902
37.1037 -7.23887
37.4546 -3.46754
37.4546 -3.76154
37.4546 -6.49579
37.4546 -7.23080
37.4797 -2.29095
37.4797 -3.05536
37.4797 -7.58303
37.5048 -2.61377
37.5048 -6.87684
37.8056 -2.60686
37.8056 -7.54613
37.8306 -2.28287
37.8306 -3.04729
37.8306 -3.42949
37.8306 -4.54671
37.8306 -4.89952
37.8306 -5.34052
37.8306 -7.19275
37.8808 -3.75174
38.1565 -7.47926
38.1816 -1.92200
38.1816 -2.27480
38.1816 -2.59821
38.1816 -3.77423
38.1816 -4.15643
38.1816 -4.53864
38.2067 -2.97984
38.2067 -3.39145
38.5325 -7.52942
38.5325 -7.94102
38.5576 -2.26616
38.5576 -3.00117
38.5576 -4.88280
38.8584 -4.87588
38.8584 -6.37530
38.8835 -1.11204
38.8835 -4.11089
38.8835 -4.49310
38.8835 -5.90432
38.8835 -7.13914
38.8835 -7.49194
38.9085 -0.78806
38.9085 -1.49367
38.9085 -3.72811
38.9085 -7.93238
38.9336 0.41793
39.1843 -6.42661
39.2094 -1.51616
39.2094 -4.48560
39.2094 -4.86781
39.2344 -1.86839
39.2344 -2.25059
39.2344 -6.01385
39.2344 -6.74886
39.2344 -7.13107
39.2595 -3.74944
39.2595 -4.13164
39.2595 -7.51270
39.2595 -7.86550
39.3096 -1.10224
39.5603 1.69656
39.5603 -6.35916
39.5854 -1.09590
39.5854 -1.50751
39.5854 -1.88972
39.5854 -2.27192
39.5854 -4.47696
39.5854 -6.74079
39.5854 -7.12300
39.5854 -7.50520
39.9363 2.61663
39.9363 -0.73503
39.9363 -1.11723
39.9363 -1.49944
39.9363 -4.08668
39.9363 -4.49829
39.9363 -4.85109
39.9363 -5.96831
39.9363 -6.35051
39.9363 -6.76212
39.9363 -7.49713
39.9614 -1.88107
39.9614 -2.26327
39.9614 -2.55728
39.9614 -7.11435
40.2873 -0.72696
40.2873 -1.10916
40.2873 -1.49137
40.2873 -1.87357
40.2873 -2.25578
40.2873 -2.54978
40.2873 -6.72465
40.2873 -7.10686
40.2873 -7.48906
40.3123 6.85894
40.3123 -3.69583
40.3123 -4.07803
40.3123 -4.46024
40.3123 -4.87184
40.3123 -5.25405
40.3123 -5.60686
40.3123 -5.98906
40.3123 -6.34187
40.6633 -0.30670
40.6633 -6.01039
40.6883 -0.71773
40.6883 -1.09994
40.6883 -2.24656
40.6883 -2.54056
40.6883 -4.86320
40.6883 -5.24540
40.6883 -5.62761
40.6883 -6.36262
40.6883 -6.74483
40.6883 -7.09763
40.6883 -7.47984
40.7134 6.86816
40.7134 6.48596
40.7134 6.10375
40.7134 -1.48157
40.7134 -1.86377
40.7134 -2.92219
40.7134 -3.33380
40.7134 -3.68660
40.7134 -4.09821
40.7385 -4.45044
41.0393 -0.29806
41.0393 -4.09071
41.0393 -4.44352
41.0393 -5.23733
41.0393 -6.00174
41.0393 -6.35455
41.0393 -6.70736
41.0393 -7.08956
41.0393 -7.47177
41.0643 6.46463
41.0643 0.08473
41.0643 -1.06189
41.0643 -1.85570
41.0643 -2.23791
41.0643 -2.56131
41.0643 -2.94352
41.0643 -3.35513
41.0643 -3.70793
41.0643 -4.85455
41.0643 -5.61896
41.0894 6.84741
41.0894 -1.47292
41.3651 0.09164
41.3902 0.47443
41.3902 -1.05440
41.3902 -4.43545
41.3902 -5.22926
41.3902 -5.99367
41.3902 -6.34648
41.3902 -6.72869
41.4153 9.91255
41.4153 0.85721
41.4153 -1.43603
41.4153 -2.55324
41.4153 -2.90605
41.4153 -3.67046
41.4153 -4.05267
41.4153 -4.84648
41.4153 -5.61089
41.4153 -7.08091
41.4153 -7.46312
41.4404 -3.34648
41.7412 9.50844
41.7412 9.12623
41.7412 1.24691
41.7412 0.89411
41.7412 0.48250
41.7412 -2.22234
41.7412 -2.54575
41.7412 -2.89855
41.7412 -3.31016
;
run;
 
DATA ALL;
SET world ob_anno;
run;
 
proc gproject DATA=ALL out=all_projected
latmin=30 latmax=55 longmin=-15 longmax=10 project=hammer dupok;
id id;
run;
 
DATA world_projected ob_anno_projected; SET all_projected;
IF flag=1 then output ob_anno_projected;
else output world_projected;
run ;
 
/*
Add the annotate variables after the grpoject (so the gproject will
be lean & fast).
*/
DATA ob_anno_projected; SET ob_anno_projected;
length FUNCTION color style $8 html $500;
xsys='2'; ysys='2'; when='a';
FUNCTION='pie'; size=0.2; rotate=360;
color='red'; style='psolid'; output;
html='title='||quote('lat: '||trim(LEFT(lat))||', '||'0d'x||'long:'||trim(LEFT(long)))||
' href='||'http://maps.google.com/maps?ie=UTF8&ll='||trim(LEFT(lat))||','||trim(LEFT(long))||'&z=8';
color='gray'; style='pempty'; output;
run;
 
%let name=map;
filename odsout '.';
goptions device=png xpixels=1200 ypixels=1200;
goptions gunit=pct ftitle="arial/bo" ftext="arial" htitle=5 htext=3 ;
 
ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" style=minimal;
 
pattern1 c=cx7f00ff v=solid;
 
proc gmap DATA=world_projected map=world_projected;
id id;
choro id / nolegend levels=1 coutline=graycc anno=ob_anno_projected
des='' name="&name";
run;
quit;