WCS用例

作者

诺尔曼巴克

联系

ittvis.com上的nbarker

作者

盖尔米林

联系

ittvis.com上的nbarker

最后更新

2019-11-21

本文档解释了如何使用MapServer通过MapServer WCS接口传递Landsat、SPOT、DEM和NetCDF临时/带状数据。感谢Steve Lime和Frank Warmerdam对这些项目的帮助。

它还演示了如何为GRIB2输出格式配置MapServer。

陆地卫星

要通过MapServer网络覆盖服务提供陆地卫星图像,请指定 OUTPUTFORMAT 对象。要获得格式支持,请安装gdal库,从命令提示符和cd到安装gdal的位置,并使用命令gdalinfo--formats。将显示所有支持格式的列表,并指定该格式是只读的**<ro>>还是读写的**<rw>>对于WCS,该格式需要支持读写(除了 GDAL's 但是,拥有wcs格式)。

在下面的例子中,陆地卫星715米分辨率的马赛克是增强压缩小波格式(ECW)。通过运行gdalinfo.exe程序,我可以验证ECW格式是否具有写权限,因此可以在mapfile中指定该格式,并使用getcoverge请求进行请求。

OUTPUTFORMAT
    NAME "ECW"
    DRIVER "GDAL/ECW"
    MIMETYPE "image/ecw"
    IMAGEMODE "BYTE"
    EXTENSION "ecw"
END
LAYER
    NAME "Landsat7"
    STATUS OFF
    TYPE RASTER
    PROCESSING "SCALE=AUTO"
    UNITS Meters
    TILEINDEX "MapServer/wcs/landsat7/l7mosaic15m.shp"
    TILEITEM "Location"
    METADATA
      "wcs_description" "Landsat 7 15m resolution mosaic"
      "wcs_name" "Landsat7"
      "wcs_label" "Landsat 7 15m resolution mosaic"
      "ows_srs" "EPSG:27700"
      "ows_extent" "0 0 700005 1050000"
      "wcs_resolution" "75 75"
      "wcs_bandcount" "3"
      "wcs_formats" "ECW"
      "wcs_enable_request" "*"
    END
 END

然后,可以通过以下元素创建URL来请求getcoverge请求(使用mapfile中设置的参数):您的服务器、mapserver程序、mapfile的位置、服务类型(wcs)、请求(getcoverge)、覆盖范围(landsat7)、bbox(0,07000051050000)、crs(epsg:27700)、resx(75)、resy(75)、format(ecw)。

SPOT

现场图像可以通过MapServer网络覆盖服务发送,类似于上面的陆地卫星示例。主要的区别在于,由于斑点是灰度图像,所以wcs波段计数为1,而不是由3个波段组成的陆地卫星图像。对于本例,将使用众所周知的geotiff格式来演示在 Mapfile 中为点数据指定什么。

OUTPUTFORMAT
    NAME "GEOTIFF"
    DRIVER "GDAL/GTiff"
    MIMETYPE "image/tiff"
    IMAGEMODE "BYTE"
    EXTENSION "tif"
END
LAYER
    NAME "SPOT"
    STATUS OFF
    TYPE RASTER
    PROCESSING "SCALE=AUTO"
    UNITS Meters
    TILEINDEX "MapServer/wcs/orthospot/spot.shp"
    TILEITEM "Location"
    METADATA
      "wcs_description" "Orthospot mosaic"
      "wcs_name" "SPOT"
      "wcs_label" "Orthospot mosaic"
      "ows_srs" "EPSG:27700"
      "ows_extent" "375960 64480 497410 200590"
      "wcs_resolution" "100 100"
      "wcs_bandcount" "1"
      "wcs_formats" "GEOTIFF"
      "wcs_nativeformat" "8-bit GeoTIF"
      "wcs_enable_request" "*"
    END
END

要在WCS映射文件中为任何数据层和格式指定的关键参数为:

- Layer Name = Create a short name for the data
- Layer Type = Raster

以下示例进一步演示了如何实现WCS,以及如何创建包含具有时间维度的层的WCS(请参见netcdf示例)。

DEM

可以通过MapServer Web覆盖服务提供16位DEM数据。

首先需要在 Mapfile 中指定输出格式

OUTPUTFORMAT
  NAME "GEOTIFFINT16"
  DRIVER "GDAL/GTiff"
  MIMETYPE "image/tiff"
  IMAGEMODE "INT16"
  EXTENSION "tif"
END

以及相应的层

LAYER
  NAME "srtm"
  STATUS OFF
  TYPE RASTER
  DATA "srtm.tif"
  PROJECTION
    "init=epsg:4326"
  END
  METADATA
    "wcs_label" "SRTM WCS TIF Server"
    "ows_extent" "-180 -90 180 90"
    "wcs_resolution" "0.00083 -0.00083"
    "ows_srs" "EPSG:4326"
    "wcs_formats" "GEOTIFFINT16"
    "wcs_nativeformat" "geotiff"
    "wcs_enable_request" "*"
  END
END

通过使用https://gdal.org/programs/gdaladdo.html中描述的gdaladdo实用程序可以提高性能

NETCDF

首先,GDAL不支持netCDF的所有版本(有很多,它是一种通用格式),因此为了保持稳定性,可能需要先将文件转换为GeoTiff格式。这可以使用netCDF库在这里实现https://www.unidata.ucar.edu/software/netCDF/。Denis Nadeau和Frank wartemdam在GDAL中添加了netCDF-CF作为只读格式,因此现在可以直接从磁盘读取CF约定的netCDF文件。

我们将z级别放置在gdal数据文件(geotiff或netcdf)的带区中,并为时间级别创建了一个形状索引。gdal数据是二维格式(x,y)和波段。netcdf是一种n-d文件格式,支持时间、x、y、z和实验参数。通过使用一组gdal netcdf/geotiff文件,可以表示这一点,并将z级别(高度)存储为数据文件中的带区。虽然是一个黑客,但是定制客户机可以通过在返回的axes描述标签中对其进行编码,从wcs的describeCoverage操作中接收关于geotiff的z-levelA带表示的重要元数据。

为了创建时间维度的形状文件,我们必须对Java代码进行一些黑客攻击,但我们也得到了它与Steve Lime的Perl脚本在MODIS MMASPerver演示下载中的工作(现在似乎不可用)。

SteveLime在modis演示中使用的Perl脚本如下,我在下面放置了内联注释。脚本假定gdaltindex已经在此目录中运行,以创建平铺索引形状和dbf文件。它假定数据文件的文件名在文件名中有日期,例如myfileyyyymmddhh.tif

 1  #!/usr/bin/perl
 2  use XBase;
 3  opendir(DIR, '.'); # open the current directory
 4  foreach $file (readdir(DIR)) {
 5    next if !($file =~ /\.dbf$/); # read the dbf file in this directory created by gdaltindex
 6    print "Working on $file...\n";
 7    $tfile = 'temporary.dbf';
 8    system("mv $file $tfile");
 9    $oldtable = new XBase $tfile or die XBase->errstr;
10    print join("\t", $oldtable->field_names) ."\n";
11    print join("\t", $oldtable->field_types) ."\n";
12    print join("\t", $oldtable->field_lengths) ."\n";
13    print join("\t", $oldtable->field_decimals) ."\n";
14    $newtable = XBase->create("name" => $file,
15                  "field_names" => [$oldtable->field_names, "IMGDATE"],  # this is the FILTERITEM in the corresponding tile index layer within your mapfile
16                  "field_types" => [$oldtable->field_types, "C"],  # character column type
17                  "field_lengths" => [$oldtable->field_lengths, 13], # length of the date string
18                  "field_decimals" => [$oldtable->field_decimals, undef]) or die "Error creating new table.";
19    foreach (0 .. $oldtable->last_record) {
20    ($deleted, @data) = $oldtable->get_record($_);
21    print "  ...record $data[0]\n";
22        # extract the date
23        $year = substr $data[0], 8, 4;  # year  is at position 8 in the filename string
24        $month = substr $data[0], 12, 2;  # month is at position 12 in the filename string
25        $day = substr $data[0], 14, 2; # day is at position 14 in the filename string
26        $hour = substr $data[0], 16, 2;  # hour is at position 16 in the filename string
27        $date =  "$year-$month-$day" . "T" . "$hour\n"; # format is YYYY-MM-DDTHH, or any ISO format
28     print "$date";
29        push @data, "$date";
30    $newtable->set_record($_, @data);
31    }
32    $newtable->close();
33    $oldtable->close();
34    unlink($tfile);
35  }

如果已经使用了Perl脚本,那么跳到下面的层定义,如果您希望自己编写代码,那么说明就在这里。

dbf文件必须具有列“location”,该列指示数据文件的位置(绝对路径或相对于映射文件位置的相对路径),以及第二列,该列可以根据需要调用,但索引时间。在我们的例子中,我们称之为“时间”:-)

然后,相应的shapefile必须包含每次带有tif文件边界框的多边形。所以ogrinfo timeindex.shp看起来像:

OGRFeature(timeIndex):116
  location(String) = mytime.tif
  time(String) = 2001-01-31T18:00:00
  POLYGON ((xxx,xxxx,.......))

将输出格式定义为

OUTPUTFORMAT
  NAME "GEOTIFF_FLOAT"
  DRIVER 'GDAL/GTiff'
  MIMETYPE 'image/tiff'
  IMAGEMODE FLOAT32
  EXTENSION 'tif'
END

然后需要在 Mapfile 中定义图块索引

LAYER
  NAME 'time_idx'
  TYPE TILEINDEX
  DATA 'timeIndex'
  FILTERITEM 'time'
  FILTER '%time%'
END

以及实际层

LAYER
  NAME 'TempData'
  STATUS OFF
  TYPE RASTER
  TILEINDEX 'time_idx'
  PROJECTION
    "init=epsg:4326"
  END
  METADATA
    "wcs_label" 'Temperature data'
    "ows_extent" '-180 -90 180 90'
    "wcs_resolution" '1.125 -1.125'
    "ows_srs" 'EPSG:4326'
    "wcs_formats" 'GEOTIFF_FLOAT'
    "wcs_nativeformat" 'netCDF'
    "wcs_bandcount" '27'
    "wcs_rangeset_axes" 'bands'
    "wcs_rangeset_label" 'Pressure (hPa units) Levels'
    "wcs_rangeset_name" 'bands'
    "wcs_rangeset_description" 'Z levels '
    "wcs_timeposition" '2001-01-01T06:00:00,2001-01-01T12:00:00,2001-01-01T18:00:00,2001-01-02T00:00:00'
    "wcs_timeitem" 'time'
    "wcs_enable_request" "*"
  END
END

tempdata覆盖层现在将允许您使用&bands=进行子集。时间=…子集参数!

要进行协调再投影,请在请求和响应中指定“CRS=ESPG:XXXX”

当您开始使用WCS和MapServer进行临时子集设置时,您会发现需要一种自动生成映射文件的方法,例如使用XSL样式表!

对于平铺索引层,您需要提供以下额外的元数据才能将其用于WCS:

"OWS_EXTENT" "10050 299950 280050 619650"
"WCS_RESOLUTION" "100 100"
"WCS_SIZE" "2700 3197"
"WCS_BANDCOUNT" "3"

如果您的图像有一个颜色表并且只有一个波段,除非您将图像模式设置为PC256而不是字节,否则它将显示灰度。

GRIB2输出格式

这要求mapserver 7.2.0用于特定格式和特定层的创建选项机制,以及gdal 2.3.0用于grib2输出支持。

将输出格式定义为

OUTPUTFORMAT
    NAME GRIB
    DRIVER "GDAL/GRIB"
    MIMETYPE "application/x-grib2"
    IMAGEMODE Float32
    EXTENSION "grib2"
    FORMATOPTION "DATA_ENCODING=SIMPLE_PACKING"
END

(查阅) GDAL GRIB driver documentation 有关可以在formatOption中提供的更多选项)

以及实际层

LAYER
    NAME temperatures
    TYPE raster
    STATUS ON
    DUMP TRUE
    DATA "temperatures.tif"

    PROJECTION
        "init=epsg:4326"
    END

    METADATA
        "ows_extent" "-180 -90 180 90"
        "wcs_label" "Test label"
        "ows_srs"   "EPSG:4326"
        "wcs_resolution" "2.4 2.4"
        "wcs_bandcount" "1"
        "wcs_imagemode" "Float32"
        "wcs_formats" "GEOTIFF_FLOAT32 GRIB"
        "wcs_description" "Test description"
        "wcs_metadatalink_href" "http://www.gdal.org/metadata_test_link.html"
        "wcs_keywordlist" "test,mapserver"
        "wcs_abstract" "abstract"

        # GDAL creation options for the GRIB output format
        # In the case where the input DATAset would be a GRIB2 product,
        # those wcs_outputformat_GRIB_creationoption_* keywords are not
        # needed, as being directly taken from the input dataset
        "wcs_outputformat_GRIB_creationoption_DISCIPLINE" "0(Meteorological)"
        "wcs_outputformat_GRIB_creationoption_IDS" "Normally this will never be used given the BAND_1_IDS override"
        "wcs_outputformat_GRIB_creationoption_BAND_1_IDS" "CENTER=54(Montreal) SUBCENTER=0 MASTER_TABLE=4 LOCAL_TABLE=0 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2017-01-05T00:00:00Z PROD_STATUS=0(Operational) TYPE=2(Analysis_and_forecast)"
        "wcs_outputformat_GRIB_creationoption_PDS_PDTN" "0"
        "wcs_outputformat_GRIB_creationoption_PDS_TEMPLATE_ASSEMBLED_VALUES" "0 0 2 47 47 0 0 1 0 103 0 2 255 -127 -2147483647"

        # WCS 2.0 metadata items
        "wcs_native_format" "application/x-grib2"
        "wcs_band_names" "Temperature_2m"
        "Temperature_2m_band_interpretation" "interp"
        "Temperature_2m_band_uom" "Celcius"
        "Temperature_2m_band_definition"       "Temperature at 2m above ground"
        "Temperature_2m_band_description"      "Temperature at 2m above ground"
        "Temperature_2m_interval" "-100 100"
        "Temperature_2m_significant_figures" "1"

        # WCS 1.x metadata items
        "wcs_nativeformat" "GRIB"
        "wcs_rangeset_axes" "bands"
        "wcs_rangeset_name" "Temperatures"
        "wcs_rangeset_label" "Bands"
        "wcs_rangeset_description" "Temperatures"
    END
END

(查阅) GDAL GRIB driver documentation “wcs_u outputformat_grib_creationoption optionname”中提供的更多选项)

学分:加拿大气象局为每层创建选项提供了资助。