Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум программистов > Fortran > Извлечение данных


Автор: spea 27.10.2012, 13:36
 Большая просьба к знающим этот язык.

 Имеется исходный файл, содержащий данные. Но в файле совершенно непонятная кодировка. AkelPad сообщает, что файл бинарный.
И есть http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/Anomaly_DOV/read_1file_min5.0 на Fortran для извлечения данных из этого файла. 

 Языки программирования почти не знаю , и разобраться как работает эта программа совершенно не получается. Пробовал её скомпилировать в Force 2.0, но выскакивает масса ошибок. По идее, в сам код нужно вносить какие-то параметры для извлечения данных из файла.

 Если вас не затруднит, посмотрите пожалуйста, что именно нужно сделать, чтобы извлечь данные.

 Исходные данные и программа находятся по адресу http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/anomalies_dov.html

 Там где написано "5 x 5-Minute Area-Mean Free-Air Gravity Anomaly Grid (Dg01_mean5.0x5.0_EGM08_to2190_WGS84_ell_nh - 36 MB)" - это исходные данные.
 "Program to Read 5-Minute Gravity Anomaly Grid (read_1file_min5.0 - 6 KB)" - текст программы.
 "Output from the Program Above (read_1file_min5.0.out01 - 7 KB)" - пример работы этой программы

Автор: Фантом 27.10.2012, 17:21
Во-первых, то, что на сайте называется "программой", на самом деле представляет собой программу на Фортране, встроенную в скрипт ksh, причем внешняя "обвязочная" часть написана, да простят меня модераторы, через задницу. Неудивительно, что скомпилировать эту химеру не удалось. Похоже, что кто-то когда-то написал код на Фортране, а уже затем кто-то другой, с Фортраном не знакомый, приделал снаружи костыли вместо переделки собственно программы.

Во-вторых, собственно фортрановский исходник также писался как минимум двумя разными людьми, причем первый из них использовал стандарт FORTRAN 77 (и делал все аккуратно), а второй - FORTRAN 90 (и опять-таки через задницу).

Соответственно, если все это зачистить, собрать и запустить, то получается следующее:
Код



               Statistics of Gravity Anomaly Values  (mGal)

       89.9583358765   0.0416666679        0.00000
       89.9583358765   0.1250000000       -0.00000
       89.9583358765   0.2083333433***************
       89.9583358765   0.2916666865        0.00000
       89.9583358765   0.3750000000       -0.00000
       89.9583358765   0.4583333433***************
       89.9583358765   0.5416666865***************
       89.9583358765   0.6250000000       -0.00000
       89.9583358765   0.7083333731        0.80194
       89.9583358765   0.7916666865       -0.00000
       89.9583358765   0.8750000000  -178477.00000
       89.9583358765   0.9583333731       -0.00489
       89.9583358765   1.0416667461***************
       89.9583358765   1.1250000000***************
       89.9583358765   1.2083333731526977024.00000
       89.9583358765   1.2916667461***************
       89.9583358765   1.3750000000120732160.00000
       89.9583358765   1.4583333731***************
       89.9583358765   1.5416667461***************

         Number of Values         20549
       Percentage of Area            -0.000
            Minimum Value   ***************
      Latitude of Minimum            -0.375
     Longitude of Minimum           263.958
            Maximum Value   ***************
      Latitude of Maximum            -0.208
     Longitude of Maximum           318.875
          Arithmetic Mean   ***************
       Area-Weighted Mean   ***************
           Arithmetic RMS   ***************
        Area-Weighted RMS   ***************
        Arithmetic S.Dev.   ***************
     Area-Weighted S.Dev.   ***************
 


Звездочки - там, где результаты не влезли в заданный формат. В принципе, если надо, я могу исправить и это...

Автор: spea 27.10.2012, 17:32
Мне на одном форуме также подсказали, что это не листинг.
Цитата
Это не фортрановский листинг, его бессмысленно так собирать. Прога под линухой. Шеловские команды формируют текстовый файл, экзешник, входной файл, все запускают и потом за собой чистят.


Цитата

Звездочки - там, где результаты не влезли в заданный формат. В принципе, если надо, я могу исправить и это...
 Сейчас тоже дошёл до этих звездочек, только мне кажется, что в программе какая-то ошибка, ибо на сайте в примере указано следующее


Код
Statistics of Gravity Anomaly Values (mGal)
 89.9583333333 0.0416666667 3.25491
 89.9583333333 0.1250000000 3.24581
 89.9583333333 0.2083333333 3.23676
 89.9583333333 0.2916666667 3.22776
 89.9583333333 0.3750000000 3.21880
 89.9583333333 0.4583333333 3.20989
 89.9583333333 0.5416666667 3.20103
 89.9583333333 0.6250000000 3.19221
 89.9583333333 0.7083333333 3.18345
 89.9583333333 0.7916666667 3.17473
 89.9583333333 0.8750000000 3.16606
 89.9583333333 0.9583333333 3.15744
 89.9583333333 1.0416666667 3.14886
 89.9583333333 1.1250000000 3.14034
 89.9583333333 1.2083333333 3.13186
 89.9583333333 1.2916666667 3.12343
 89.9583333333 1.3750000000 3.11505
 89.9583333333 1.4583333333 3.10672
 89.9583333333 1.5416666667 3.09843

 Number of Values 9331200
 Percentage of Area 100.000
 Minimum Value -361.775
 Latitude of Minimum 19.375
 Longitude of Minimum 293.542
 Maximum Value 869.094
 Latitude of Maximum 10.792
 Longitude of Maximum 286.292
 Arithmetic Mean -0.471
 Area-Weighted Mean 0.000
 Arithmetic RMS 32.899
 Area-Weighted RMS 34.457
 Arithmetic S.Dev. 32.896
 Area-Weighted S.Dev. 34.457

Автор: Фантом 27.10.2012, 18:35
Скорее не ошибка, а особенность интерпретации двоичных данных. Собственно, судя по тексту программы, она должна считывать из файла с данными 2160*4320 вещественных чисел, по 4 байта на число. Проблема в том, что на разных платформах форма хранения вещественных чисел может быть разной (у меня в итоге получаются отсчеты, различающиеся порядков на 60, что, наверное, не слишком осмысленно). 

Можно попробовать подобрать форму, однако для этого нужно знать, какими примерно по величине могут быть данные. Вы можете это сообщить?

Автор: spea 27.10.2012, 18:47
Цитата

Можно попробовать подобрать форму, однако для этого нужно знать, какими примерно по величине могут быть данные. Вы можете это сообщить?
 Не совсем понял вопрос. Значения первого столбика от 90 до -90. Значения второго от 0 до 360. Для каждого шаг 0,083333333. В третьем столбике значения не подчиняются какому-либо закону и варьируются где-то от -100 до +100 без резких перепадов.
Или общее количество строк данных нужно знать? 3110400 строк, что соответствует 2160*4320 разделить на 3 столбика.

Если скомпилировать этот код, то получается почти как в примере, только в третьей колонке сначала 0,00000, затем данные смещены вниз на одну строчку.
Код
      program readmin1
      implicit real*8(a-h,o-z)
      parameter(exclud=9999.d0,
     $   dlatg=5.0d0/60.d0,dlong=5.0d0/60.d0,nrowsg=02160,ncolsg=04320)
         real*4 data(ncolsg)
      dimension sdat(ncolsg,2),stati(22)
c-----------------------------------------------------------------------
c
c   Read-in Gravity Anomaly 5.0x5.0 minute mean value file.
c
open (unit   = 1,
     &     file   = 'fort.1',
     &     access = 'stream',
 &     convert='big_endian')
      write(6,6001)
 6001 format(///15x,'Statistics of Gravity Anomaly Values  (mGal)',/)
      do i = 1, nrowsg
        read(1) (  data(j),j=1,ncolsg)
        rlat = 90.d0 - (i-1.d0)*dlatg - dlatg*0.5d0
        do j = 1, ncolsg
          clon = (j-1.d0)*dlong + dlong*0.5d0
c
c   Print a few values.
c
          if(i.eq.   1.and.j.lt.20) then  !  top row
            write(6,6101) rlat,clon,data(j)
 6101       format(5x,2f15.10,f15.5)
          endif  !  top row
          sdat(j,1) = data(j)
          sdat(j,2) = 1.d0
        enddo  !  j
        call stats(90.d0,0.d0,i,dlatg,dlong,nrowsg,ncolsg,sdat,exclud,0,
     $             stati)
      enddo  !  i
c-----------------------------------------------------------------------
      stop
      end
C
C
C
      SUBROUTINE STATS(TOPLAT,WSTLON,I,GRDN,GRDE,NROWS,NCOLS,DATA,
     $                 EXCLUD,ISIG,STAT)
C-----------------------------------------------------------------------
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*20 DLABEL(14)
      CHARACTER*20 SLABEL( 8)
      DIMENSION DATA(NCOLS,2)
      DOUBLE PRECISION DEXCLUD,DG,SD,STAT(22)
      SAVE
      DATA PI/3.14159265358979323846D+00/
      DATA DLABEL/'    Number of Values','  Percentage of Area',
     $            '       Minimum Value',' Latitude of Minimum',
     $            'Longitude of Minimum','       Maximum Value',
     $            ' Latitude of Maximum','Longitude of Maximum',
     $            '     Arithmetic Mean','  Area-Weighted Mean',
     $            '      Arithmetic RMS','   Area-Weighted RMS',
     $            '   Arithmetic S.Dev.','Area-Weighted S.Dev.'/
      DATA SLABEL/'       Minimum Sigma',' Latitude of Minimum',
     $            'Longitude of Minimum','       Maximum Sigma',
     $            ' Latitude of Maximum','Longitude of Maximum',
     $            'Arithmetic RMS Sigma','Area-wghtd RMS Sigma'/
C-----------------------------------------------------------------------
      IF(I.EQ.1) THEN
      DEXCLUD=EXCLUD
      DTR=PI/180.D0
      FOURPI=4.D0*PI
      DPR=GRDN*DTR
      DLR=GRDE*DTR
      CAREA=2.D0*DLR*SIN(DPR/2.D0)
      DO 10 K=1,22
      STAT(K)=0.D0
   10 CONTINUE
      STAT( 3)= DEXCLUD
      STAT( 6)=-DEXCLUD
      STAT(15)= DEXCLUD
      STAT(18)= 0.D0
      ENDIF
C-----------------------------------------------------------------------
      DLAT=TOPLAT-(I-1.D0)*GRDN-GRDN/2.D0
      COLATC=(90.D0-DLAT)*DTR
      AREA=CAREA*SIN(COLATC)
C-----------------------------------------------------------------------
      DO 110 J=1,NCOLS
      DLON=WSTLON+(J-1.D0)*GRDE+GRDE/2.D0
      DG=DATA(J,1)
      SD=DATA(J,2)
      IF(DG.LT.DEXCLUD) THEN
C-----------------------------------------------------------------------
      STAT( 1)=STAT( 1)+1.D0
      STAT( 2)=STAT( 2)+AREA
      IF(DG.LE.STAT( 3)) THEN
      STAT( 3)=DG
      STAT( 4)=DLAT
      STAT( 5)=DLON
      ENDIF
      IF(DG.GE.STAT( 6)) THEN
      STAT( 6)=DG
      STAT( 7)=DLAT
      STAT( 8)=DLON
      ENDIF
      STAT( 9)=STAT( 9)+DG
      STAT(10)=STAT(10)+DG*AREA
      STAT(11)=STAT(11)+DG**2
      STAT(12)=STAT(12)+DG**2*AREA
      IF(SD.LE.STAT(15)) THEN
      STAT(15)=SD
      STAT(16)=DLAT
      STAT(17)=DLON
      ENDIF
      IF(SD.GE.STAT(18)) THEN
      STAT(18)=SD
      STAT(19)=DLAT
      STAT(20)=DLON
      ENDIF
      STAT(21)=STAT(21)+SD**2
      STAT(22)=STAT(22)+SD**2*AREA
C-----------------------------------------------------------------------
      ENDIF
  110 CONTINUE
C-----------------------------------------------------------------------
      IF(I.NE.NROWS) RETURN
      IF(STAT(1).GT.0.D0) THEN
      STAT( 9)=STAT( 9)/STAT( 1)
      STAT(10)=STAT(10)/STAT( 2)
      STAT(11)=SQRT(STAT(11)/STAT( 1))
      STAT(12)=SQRT(STAT(12)/STAT( 2))
      STAT(13)=SQRT(STAT(11)**2-STAT( 9)**2)
      STAT(14)=SQRT(STAT(12)**2-STAT(10)**2)
      STAT(21)=SQRT(STAT(21)/STAT( 1))
      STAT(22)=SQRT(STAT(22)/STAT( 2))
      STAT( 2)=STAT( 2)/FOURPI*100.D0
      ELSE
      DO 120 J=3,22
      STAT(J)=DEXCLUD
  120 CONTINUE
      ENDIF
C=======================================================================
      NUM=INT(STAT(1))
      WRITE(6,6001) DLABEL(1),NUM
 6001 FORMAT(/5X,A20,3X,I11)
      DO 210 K=2,14
      WRITE(6,6002) DLABEL(K),STAT(K)
 6002 FORMAT(5X,A20,3X,F15.3)
  210 CONTINUE
      WRITE(6,6003)
 6003 FORMAT(' ')
      IF(ISIG.EQ.1) THEN
      DO 220 K=1,8
      WRITE(6,6002) SLABEL(K),STAT(K+14)
  220 CONTINUE
      ENDIF
C=======================================================================
      RETURN
      END

Автор: Фантом 27.10.2012, 19:39
Цитата(spea @  27.10.2012,  19:47 Найти цитируемый пост)
 Не совсем понял вопрос. Значения первого столбика от 90 до -90. Значения второго от 0 до 360. Для каждого шаг 0,083333333. В третьем столбике значения не подчиняются какому-либо закону и варьируются где-то от -100 до +100 без резких перепадов.


Я не про вывод. Я про значения во входном файле. Какими они могут быть?

Автор: spea 27.10.2012, 19:48
как
Цитата(Фантом @  27.10.2012,  19:39 Найти цитируемый пост)
Я не про вывод. Я про значения во входном файле. Какими они могут быть?
Наверно я дико туплю, но в файле должны быть данные такие же как и в примере

 89.9583333333  0.0416666667  3.25491
 89.9583333333  0.1250000000  3.24581
 89.9583333333  0.2083333333  3.23676
 89.9583333333  0.2916666667  3.22776
 89.9583333333  0.3750000000  3.21880
 89.9583333333  0.4583333333  3.20989
 89.9583333333  0.5416666667  3.20103
 89.9583333333  0.6250000000  3.19221
 89.9583333333  0.7083333333  3.18345
 89.9583333333  0.7916666667  3.17473

Автор: spea 27.10.2012, 20:49
Вычислил, что в начале должно быть так
Код
      open (unit   = 1,
     &     file   = 'fort.1',
     &     form   = 'unformatted',
     &     access = 'sequential',
     &     convert='big_endian')



Открывается также как и в примере. Вопрос в том, как записать данные в нормальном формате в текстовый файл?

Автор: Фантом 27.10.2012, 21:07
Цитата(spea @  27.10.2012,  20:48 Найти цитируемый пост)
Наверно я дико туплю, но в файле должны быть данные такие же как и в примере


Точнее, в последнем столбце примера. Т.е. что-то порядка единицы?

Тогда байты надо расставлять в обратном порядке.  Получается то, что в примере, но опять-таки со сдвигом 3-го столбца на одну позицию. Если считать, что первые четыре байта в данных пустые (на что очень похоже - файл чуть длиннее, чем требуется), то получается соответствие.
Вот простенькая программа на C99, которая перелопачивает исходный файл (переименованный в  data.dat) в текстовое представление:
Код

#include<stdio.h>

int main()
{
    const int nrowsg=02160;
    const int ncolsg=04320;
    const double dlatg=5.0/60.0;
    const double dlong=5.0/60.0;
    double rlat,clon;

    FILE *fi = fopen("data.dat","rb");
    FILE *fo = fopen("result.dat","w+");

    char by[4];
    float re;
    for(int k=3;k>-1;--k) by[k]=getc(fi);
    
    for(int i = 0;i<nrowsg;++i) for(int j=0;j<ncolsg;++j)
    {
        rlat = 90.0 - i*dlatg - dlatg*0.5;
        clon = j*dlong + dlong*0.5;
        for(int k=3;k>-1;--k) by[k]=getc(fi);
        re=*(float*)by;
        fprintf(fo,"%lf  %lf  %f\n",rlat,clon,re);
    }
    fclose(fi);
    fclose(fo);
}


Если хотите, могу и результат выслать, но напишите, куда - он и в ужатом виде больше 20М.

Автор: spea 27.10.2012, 21:19
Цитата

Вот простенькая программа на C99, которая перелопачивает исходный файл (переименованный в  data.dat) в текстовое представление:

Только что её скомпилировал. В результирующем файле почему-то не все точки. Должно быть 3110400 строк, а получается 2562816. 


Напутал я с числом строк, их как раз должно быть 9331200=2160*4320

Автор: Фантом 27.10.2012, 21:53
Цитата(spea @  27.10.2012,  22:19 Найти цитируемый пост)
В результирующем файле почему-то не все точки. Должно быть 3110400 строк, а получается 2562816.  


Откуда? Там в цикле читается 2160*4320 блоков по четыре байта. Соответственно, должно получаться 9331200 строк (столько, собственно, и получается).

Автор: spea 27.10.2012, 21:54
Немного подправил код
Код
#include<stdio.h>

int main()
{
    const int nrowsg=2160;
    const int ncolsg=4320;
    const double lath=5.0/60.0;
    const double lonh=5.0/60.0;
    double lat, lon;

    FILE *fi = fopen("d:\\data.dat","rb");
    FILE *fo = fopen("d:\\result.txt","w+");

    char by[4];
    float faa;
    for(int k=3;k>-1;--k) by[k]=getc(fi);
    
    for(int i = 0;i<nrowsg;++i)
    { 
    for(int j=0;j<ncolsg;++j)
    {
        lon = 90.0 - i*lonh - lonh*0.5;
        lat = j*lath + lath*0.5;
        for(int k=3;k>-1;--k) by[k]=getc(fi);
        faa=*(float*)by;
        fprintf(fo,"%lf  %lf  %f\n",lon,lat,faa);
    }
    }
    fclose(fi);
    fclose(fo);
}




Теперь строчек, сколько нужно. 

Фантом, огромнейшее спасибо! Вы меня невероятно выручили.

Добавлено @ 22:01
Цитата

Если хотите, могу и результат выслать, но напишите, куда - он и в ужатом виде больше 20М.

Вот только у меня получается заархивированный 50M, не сжатый 300M. 

Автор: Фантом 27.10.2012, 22:07
Цитата(spea @  27.10.2012,  22:54 Найти цитируемый пост)
Немного подправил код

По идее разницы быть не должно. Может быть, у Вас оптимизация компилятора хулиганит?

Цитата(spea @  27.10.2012,  22:54 Найти цитируемый пост)

Вот только у меня получается заархивированный 50M, не сжатый 300M.  

Аналогично. Я просто в первый раз прервал программу раньше и ужимал только то, что успело получиться к тому моменту. Потому и написал "больше".

Автор: spea 27.10.2012, 22:12
Цитата

Может быть, у Вас оптимизация компилятора хулиганит?
 Может быть, я в этом не силён. Пользуюсь Dev-C++.

Автор: spea 28.10.2012, 00:58
Столкнулся с такой проблемой. 
Через каждые 4320 точек встречается 2 пары нулей в третьем столбце. И остальные данные получаются смещёнными на 2 строчке и это смещение накапливается.
Вот пример резкого перепада:
Код
89.791667  359.791667  11.624873
89.791667  359.875000  11.579295
89.791667  359.958333  11.533430
89.708333  0.041667  11.487279
89.708333  0.125000  11.440846
89.708333  0.208333  11.394135
89.708333  0.291667  11.347147
89.708333  0.375000  0.000000
89.708333  0.458333  0.000000
89.708333  0.541667  16.022243
89.708333  0.625000  15.987768
89.708333  0.708333  15.952506



Цитата
Если считать, что первые четыре байта в данных пустые (на что очень похоже - файл чуть длиннее, чем требуется), то получается соответствие.
 Да, размер исходного файла 2160*4320*4 + 4320*4. Каждые 4322 чисел начинаются и заканчиваются нулём: 
0.00000
далее 4320 нужных значений
0.00000
0.00000
далее 4320 нужных значений
0.00000

и т.д.

Автор: spea 28.10.2012, 01:56
Фантом, проверьте пожалуйста, получается ли у вас без нулей?

upd. Немного коряво, но итогового результата добился. Считал исходный файл "как есть" на 2160*4320+4320 строк, не создавая первые два столбца. Затем удалил 4320 нулей. Занёс полученные данные в массив и присоединил его к первым двум столбцам. Получилось то, что нужно: всё совпало, смещений и нулей нет. 
Буду рад, если подскажите более простое решение. Хотя, в принципе, можно считать что проблема решена.

P.S. Извиняюсь за оверпостинг. Думал, что сообщение приклеится к предыдущему.

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)