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”中提供的更多选项)
学分:加拿大气象局为每层创建选项提供了资助。