修改USGS静态地形数据方法探讨

天目神眉

先来点介绍：
摘译自USGS文件Standards for Digital Elevation Models)http://edc.usgs.gov/
30秒地形数据可以从USGS的匿名ftp站点edcftp.cr.usgs.gov下的目录：/pub/data/gtopo30/global中获得,其中e100n40为中国区.30'大约1km*1km网格.
1、概述
    USGS 是美国地质调查局（U.S. Geological Survey）的英文缩写，USGS负责管理美国全国的数字地图数据的采集与分发。
1．1 USGS DEM数据产品的种类
（1）7.5－分DEM （一般采用30米格网间距，UTM投影），水平格网间距可以去1-30之间任意整数。DEM的范围大约为无重叠的标准的USGS 7.5分地理格网。
（2）30－分DEM（2×2秒间距）DEM的范围大约为无重叠的标准的USGS 30分×30分地理格网。
（3）1－度DEM（2×2秒间距）DEM的范围大约为无重叠的标准的USGS 1度×1度地理格网。
（4）7.5－分阿拉斯加DEM（1×2秒间距，纬度，经度），范围与7.5－分DEM基本相同除了在经度元素变化从最南端的10分变化至最北端的18分。
（5）15-分阿拉斯加DEM（2×3秒间距，纬度，经度），在阿拉斯加最南端的覆盖范围为15分（纬度）×20分（经度），在最北端经度范围变为36分。
1．2 USGS DEM的格式
    USGS DEM文件由逻辑记录A、B、C组成，其中第一部分是文件头记录 type A，主要记录了DEM数据有关的信息；第二部分是断面数据type B，分为断面头数据和DEM数据实体；第三部分是精度信息type C,可以省略。
    USGS DEM数据以ASCII码形式存储，逻辑记录A、B、C格式说明分别见附表1、2、3。逻辑记录A、B、C都以1024字节长度作为逻辑记录单位，不足1024的用空格补齐。逻辑记录B通常包含多个1024字节长度的逻辑记录单位。为了有效利用空间每4个逻辑记录单位组成一个物理记录单位（4096 字节）。
2、DEM 的数据结构
    USGS DEM主要采用两种类型的格网，采用UTM投影和采用地理坐标以秒为单位的格网。这里主要介绍以秒为单位的格网数据结构。一个典型的以秒为单位的DEM数据结构如图1-2所示，数据覆盖区域是一个地理上的矩形。DEM数据的四个角点坐标记录在逻辑记录A中，详见附表1。每一个断面的起始点坐标记录在逻辑记录B中，详见附表2。这些坐标描述了DEM数据的矩形形状和每个断面的起始点坐标。以上关于秒制DEM的规定适用于除7.5分UTM DEM以外所有DEM数据。
3、USGS DEM质量控制信息
Level 1
        Level 1DEM 数据通常采用标准记录格式，数据通常是7.5－分DEM，数据通常来源于航片采集，通常DEM数据要求均方根误差RMSE不应超过7米，最大不超过15米，最大误差不超过50米。
Level 2
        Level 2DEM 数据通常经过了编绘，最大允许均方根误差为等高线间距的1/2，最大误差为1个高线间距。
Level 3
        Level 3DEM 数据通常来源于线划图，最大允许均方根误差为等高线间距的1/3，最带误差为等高线间距的2/3。

而模式中用的只有10m,5m,2m,30s 四种，看完介绍，再来看看ARWUsersGuide中对修改静态地形数据的相关说明；
用Geogrid二进制格式写静态数据

geogrid.exe 用来插值的地形静态数据以标准的2-D和3-D数值形式写成一个简单的二进制格式。当用户有了新的静态数据资源可以通过把数据写成二进制形式让WPS来使 用。geogrid格式支持单层或多层连续的场，场的分类（作为区域的分类）及为每种分类给出的一小部分的场。根据二进制格式的含义，下面给出了一个 30-second的USGS土地利用数据集及的模板：

对 于一个主导分类的分类要素场，数据首先要存储在一个规则的2-D整型数组中，每个整数都给出了与格点相对应的主导分类的说明。考虑到这个数组的作用，数据 是由低（空）到高（空）或由南向北一行行的写入文件中。例如，在上边的图中，n×m的数组的各个元素是按照x11, x12, ..., x1m, x21, ... , x2m, ... , xn1, ... , xnm的顺序写入的。尽管可以通过调整“index”文件中的参数endian=little来设置数据集中排列顺序为little-endian，在写 入文件时，每个元素是按big-endian的顺序排列的（例如，对于4-字节整数ABCD，字节A被存储在最低的地址中，而字节D被存储在最高的字节 中）。而且虽然在一个数组中用所需的最小数字来代表一个完整的数值范围是最好的，但在这里一个文件中必须用相同的为位数来存储各个元素。当写入二进制数据 时，不用头文件，记录标志，或者多余的位数。例如，一个两位的1000×1000数组文件的大小是2,000,000。因为Fortran无格式的写记录 标志，因此它是不可能直接写出一个geogrid的二进制文件的，但是可以通过用C或Fortran在写或读时调用read_geogrid.c和 write_geogrid.c（在geogrid/src目录下）来实现。与主导分类类似的是连续或真实的要素场。像主导分类一样，单层连续场作为2- D数组首先被组织，然后一行行写入二进制文件。尽管如此，因为一个连续长可以包含非整型或负值，所以每个文件中的存储表示有点复杂。数组中所有的元素都要 先被转成整型。首先通过扫描所有的元素，找出有需要的精度的，然后把剩下的部分通过四舍五入来移除。例如，要求精度达到小数点后三位，则-2.71828 需要被舍入成-2.718。然后所有的的数组元素都被转化成整数，如果有负数，则会有第二次转化：如果用一个字节储存一个元素，则每个负数都会添加28； 如果是两个字节则是216；三个字节是224；四个字节是232。这对无需转化成正数的过程很重要。最后的正的整数数组按照主导分类的格式被写入文件。多 层的联系场类似单层的。对于一个n × m × r数组，先按照上边的方法转化成正的整数场，然后从低到高开始写每个层次。分类场可以被看成是一个多层连续场的一部分，它的层次是k（1 ≤ k ≤ r）。当把一个场按照geogrid二进制格式写入一个文件时，用户应该遵守geogrid程序的命名规则，即希望数据文件的名字格式是xstart- xend.ystart-yend，其中xstart, xend, ystart,和yend分别是5位正整数，数值的起止x指针、起止y指针都包含在文件中；指针起始为1，而不是0。所以，一个800×1200的数组 （如800行，1200列）可以被命名为00001-01200.00001-00800。当一个数据被分成多个部分，每个部分都是一个规则的矩形数组， 每个数组都是一个单独的文件。在这种情况下，数组相对的位置是由每个数组的文件名中的x和y指针决定的。尽管如此，值得注意的是，数据集中的每个部分都要 有相同的x-和y-维数，并且每一块数据不能相互重叠；并且每块的起止都要在索引的范围内。例如，全球的30-second的USGS地形数据集被分成多 个1200×1200的数组，每个数组都代表了一个10°×10°的区域；这个文件的最西南角的位置是（90S, 180W），命名为00001-01200.00001-01200，而最东北角的位置是（90N, 180E），命名为42001-43200.20401-21600。如果数据集被分成多个小的部分，而且格点的数字不能在x-方向上被平分，那最后剩下 的几列数据就要用标志数（index file文件中参数missing_value确定）来填补。例如，一个x-方向上有2456个格点的数据集，如果分成3个部分的话，在规定范围内每块应 该是1–820,821–1640,和1641–2460，而2457-2460之间就要用标志数（类似缺测值）来填充。
因为起 止索引必须是5位数字，所以一个场是不能在x-或y-方向上超过99999个数据格点的，如果超过了，用户要把这个数据集分成多个geogrid能识别的 小的部分。例如，一个大的全球数据集可以把东半球和西半球分成多个小的部分。除了二进制数据文件，geogrid要求每一个数据集都有一个额外的说明文 件，即“index”文件，因此，两个数据集是不能放在同一个目录下的。其实，当要处理数据集时，index文件是最先被geogrid检索的，因为这个 文件里包含了所有可能与用到的数据相关的有用信息，下面是landuse的一个index的模板：
type=categorical
category_min=1
category_max=24
projection=regular_ll
dx=0.00833333
dy=0.00833333
known_x=1.0
known_y=1.0
known_lat=-89.99583
known_lon=-179.99583
wordsize=1
tile_x=1200
tile_y=1200
tile_z=1
units="category"
description="24-category USGS landuse"



       诸如00001-01200.00001-01200（1200×1200=1440000byte）命名的文件很多，分区域存储。可以很方便的查看read_geogrid.c和 write_geogrid.c两个程序，由于是用C写的，对于好久没用C的人来说，真的是没看明白。既然是二进制方式存储的，又附有相关的说明，试了一段时间，竟然没有成功。好吧，看来它的存储真的不一般了。
    欢迎大家做过相关工作的朋友来说说自己的看法！
    分别附上read_geogrid.c和 write_geogrid.c两个程序：
/* File: read_geogrid.c
   Sample subroutine to read an array from the geogrid binary format.
   Notes: Depending on the compiler and compiler flags, the name of
   the read_geogrid() routine may need to be adjusted with respect
   to the number of trailing underscores when calling from Fortran.
   Michael G. Duda, NCAR/MMM
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#ifdef _UNDERSCORE
#define read_geogrid read_geogrid_
#endif
#ifdef _DOUBLEUNDERSCORE
#define read_geogrid read_geogrid__
#endif
#define BIG_ENDIAN    0
#define LITTLE_ENDIAN 1
int read_geogrid(
      char * fname,            /* The name of the file to read from */
      int * len,               /* The length of the filename */
      float * rarray,          /* The array to be filled */
      int * nx,                /* x-dimension of the array */
      int * ny,                /* y-dimension of the array */
      int * nz,                /* z-dimension of the array */
      int * isigned,           /* 0=unsigned data, 1=signed data */
      int * endian,            /* 0=big endian, 1=little endian */
      float * scalefactor,     /* value to multiply array elements by before truncation to integers */
      int * wordsize,          /* number of bytes to use for each array element */
      int * status)
{
   int i, ival, cnt, narray;
   int A2, B2;
   int A3, B3, C3;
   int A4, B4, C4, D4;
   unsigned char * c;
   char local_fname[1024];
   FILE * bfile;
   *status = 0;
   narray = (*nx) * (*ny) * (*nz);
   /* Make a null-terminated local copy of the filename */
   strncpy(local_fname,fname,*len);
   local_fname[*len]='\0';
   /* Attempt to open file for reading */
   if (!(bfile = fopen(local_fname,"rb")))
   {
      *status = 1;
      return 1;
   }
   /* Allocate memory to hold bytes from file and read data */
   c = (unsigned char *)malloc(sizeof(unsigned char)*(*wordsize) * narray);
   cnt = fread((void *)c, sizeof(unsigned char), narray*(*wordsize), bfile);

   fclose(bfile);
   if (cnt == 0)
   {
      *status = 1;
      return 1;
   }
   /*
      Set up byte offsets for each wordsize depending on byte order.
      A, B, C, D give the offsets of the LSB through MSB (i.e., for
      word ABCD, A=MSB, D=LSB) in the array from the beginning of a word
   */
   if (*endian == BIG_ENDIAN) {
      A2 = 0; B2 = 1;
      A3 = 0; B3 = 1; C3 = 2;
      A4 = 0; B4 = 1; C4 = 2; D4 = 3;
   }
   else {
      B2 = 0; A2 = 1;
      C3 = 0; B3 = 1; A3 = 2;
      D4 = 0; C4 = 1; B4 = 2; A4 = 3;
   }
   /* Convert words from native byte order */
   switch(*wordsize) {
      case 1:
         for(i=0; i<narray; i++)
         {
            ival = (int)(c);      
            if ((*isigned) && (ival > (1 << 7))) ival -= (1 << 8);
            rarray = (float)ival;
         }
         break;
      case 2:
         for(i=0; i<narray; i++)
         {
            ival = (int)((c[2*i+A2]<<8) | (c[2*i+B2]));      
            if ((*isigned) && (ival > (1 << 15))) ival -= (1 << 16);
            rarray = (float)ival;
         }
         break;
      case 3:
         for(i=0; i<narray; i++)
         {
            ival = (int)((c[3*i+A3]<<16) | (c[3*i+B3]<<8) | c[3*i+C3]);      
            if ((*isigned) * (ival > (1 << 23))) ival -= (1 << 24);
            rarray = (float)ival;
         }
         break;
      case 4:
         for(i=0; i<narray; i++)
         {
            ival = (int)((c[4*i+A4]<<24) | (c[4*i+B4]<<16) | (c[4*i+C4]<<8) | c[4*i+D4]);      
            if ((*isigned) && (ival > (1 << 31))) ival -= (1 << 32);
            rarray = (float)ival;
         }
         break;
   }
   free(c);
   /* Scale real-valued array by scalefactor */
   if (*scalefactor != 1.0)
   {
      for (i=0; i<narray; i++)
         rarray = rarray * (*scalefactor);
   }
   return 0;
}

/* File: write_geogrid.c
   Sample subroutine to write an array into the geogrid binary format.
   Side effects: Upon completion, a file named 00001-<NX>.00001-<NY> is
   created, where <NX> is the argument nx and <NY> is the argument ny,
   both in i5.5 format.
  
   Notes: Depending on the compiler and compiler flags, the name of
   the write_geogrid() routine may need to be adjusted with respect
   to the number of trailing underscores when calling from Fortran.
   Michael G. Duda, NCAR/MMM
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#ifdef _UNDERSCORE
#define write_geogrid write_geogrid_
#endif
#ifdef _DOUBLEUNDERSCORE
#define write_geogrid write_geogrid__
#endif
#define BIG_ENDIAN    0
#define LITTLE_ENDIAN 1
int write_geogrid(
      float * rarray,          /* The array to be written */
      int * nx,                /* x-dimension of the array */
      int * ny,                /* y-dimension of the array */
      int * nz,                /* z-dimension of the array */
      int * isigned,           /* 0=unsigned data, 1=signed data */
      int * endian,            /* 0=big endian, 1=little endian */
      float * scalefactor,     /* value to divide array elements by before truncation to integers */
      int * wordsize )         /* number of bytes to use for each array element */
{
   int i, narray;
   int A2, B2;
   int A3, B3, C3;
   int A4, B4, C4, D4;
   unsigned int * iarray;
   unsigned char * barray;
   char fname[24];
   FILE * bfile;
   narray = (*nx) * (*ny) * (*nz);
   iarray = (unsigned int *)malloc(sizeof(int) * narray);
   barray = (unsigned char *)malloc(sizeof(unsigned char) * narray * (*wordsize));
   /* Scale real-valued array by scalefactor and convert to integers */
   for (i=0; i<narray; i++)
      iarray = (unsigned int)(rarray / (*scalefactor));
   /*
      Set up byte offsets for each wordsize depending on byte order.
      A, B, C, D give the offsets of the LSB through MSB (i.e., for
      word ABCD, A=MSB, D=LSB) in the array from the beginning of a word
   */
   if (*endian == BIG_ENDIAN) {
      A2 = 0; B2 = 1;
      A3 = 0; B3 = 1; C3 = 2;
      A4 = 0; B4 = 1; C4 = 2; D4 = 3;
   }
   else {
      B2 = 0; A2 = 1;
      C3 = 0; B3 = 1; A3 = 2;
      D4 = 0; C4 = 1; B4 = 2; A4 = 3;
   }
   /* Place words into storage byte order */
   switch(*wordsize) {
      case 1:
         for(i=0; i<narray; i++) {
            if (iarray < 0 && *isigned) iarray += (1 << 8);
            barray[(*wordsize)*i] = (unsigned char)(iarray & 0xff);
         }
         break;
      case 2:
         for(i=0; i<narray; i++) {
            if (iarray < 0 && *isigned) iarray += (1 << 16);
            barray[(*wordsize)*i+A2] = (unsigned char)((iarray >> 8) & 0xff);
            barray[(*wordsize)*i+B2] = (unsigned char)( iarray       & 0xff);
         }
         break;
      case 3:
         for(i=0; i<narray; i++) {
            if (iarray < 0 && *isigned) iarray += (1 << 24);
            barray[(*wordsize)*i+A3] = (unsigned char)((iarray >> 16) & 0xff);
            barray[(*wordsize)*i+B3] = (unsigned char)((iarray >> 8)  & 0xff);
            barray[(*wordsize)*i+C3] = (unsigned char)( iarray        & 0xff);
         }
         break;
      case 4:
         for(i=0; i<narray; i++) {
            if (iarray < 0 && *isigned) iarray += (1 << 32);
            barray[(*wordsize)*i+A4] = (unsigned char)((iarray >> 24) & 0xff);
            barray[(*wordsize)*i+B4] = (unsigned char)((iarray >> 16) & 0xff);
            barray[(*wordsize)*i+C4] = (unsigned char)((iarray >> 8)  & 0xff);
            barray[(*wordsize)*i+D4] = (unsigned char)( iarray        & 0xff);
         }
         break;
   }
   sprintf(fname,"%5.5i-%5.5i.%5.5i-%5.5i",1,*nx,1,*ny);
   /* Write array to file */
   bfile = fopen(fname,"wb");
   fwrite(barray,sizeof(unsigned char),narray*(*wordsize),bfile);
   fclose(bfile);
   free(iarray);
   free(barray);
   return 0;
}

  最后上传一个landuse的文件，方便有兴趣的朋友尝试。
00001-01200.00001-01200 (1.37 MB, 下载次数: 101)

2012-4-8 20:00 上传

点击文件名下载附件

zs4520

想请教楼主运行write_geogrid.c时出现以下提示，是什么原因。运行过程中没有出错，但最后没生成文件
write_geogrid.c:27:1: warning: "BIG_ENDIAN" redefined
In file included from /usr/include/bits/waitstatus.h:65,
                 from /usr/include/stdlib.h:43,
                 from write_geogrid.c:16:
/usr/include/endian.h:47:1: warning: this is the location of the previous definition
write_geogrid.c:28:1: warning: "LITTLE_ENDIAN" redefined
/usr/include/endian.h:46:1: warning: this is the location of the previous definition
write_geogrid.c: In function ?.rite_geogrid?.
write_geogrid.c:102: warning: left shift count >= width of type

随缘

最近想做关于地形的一些东西，天目这贴及时雨啊

冬日情愫+

，谢谢楼主哦！

xiaoaiai

哇 天目同学真牛~~好详细哇

夏天爱裤衩

这个不错，大神啊

天目神眉

随缘 发表于 2012-4-8 20:09

                               
登录/注册后可看大图

最近想做关于地形的一些东西，天目这贴及时雨啊

好像在群里讨论过哦，不过还有很多问题还没解决，共同期待

wangyaamy

我才刚改完geo,你已经改到USGS了！支持天目，真心佩服你的探索精神！也非常佩服你的慷慨大度！高度赞扬{:eb502:}{:eb502:}{:eb502:}{:eb502:}

wangyaamy

等我运行完这个个例，我也来试试替换USGS或者MODIS资料

随缘

wangyaamy 发表于 2012-4-9 10:35

                               
登录/注册后可看大图

我才刚改完geo,你已经改到USGS了！支持天目，真心佩服你的探索精神！也非常佩服你的慷慨大度！高度赞扬{:eb ...

改地形有实际意义吗？好多人在质疑改地形好假，改变地形的那段高度的气压怎么补偿？

wangyaamy

不好意思，我不是地形的，我是修改植被的特征的