關閉
標題：dem grd hdr 轉 tiff mapserver 發地形圖
內容：
開放資料 from :
https://www.tgos.tw/MDE/MetaData/CRUD/ViewMetaData?GUID=72047ba9-c161-442e-b641-904220c740b4&VIEW_TYPE=only&SOURCE=1
<?php
// grd_to_tif.php
function my_pathinfo($path)
{
$dirname = '';
$basename = '';
$extension = '';
$filename = '';
$pos = strrpos($path, '/');
if($pos !== false) {
$dirname = substr($path, 0, strrpos($path, '/'));
$basename = substr($path, strrpos($path, '/') + 1);
} else {
$basename = $path;
}
$ext = strrchr($path, '.');
if($ext !== false) {
$extension = substr($ext, 1);
}
$filename = $basename;
$pos = strrpos($basename, '.');
if($pos !== false) {
$filename = substr($basename, 0, $pos);
}
return array (
'dirname' => $dirname,
'basename' => $basename,
'extension' => $extension,
'filename' => $filename
);
}
function subname($fname){
//$pathinfo=pathinfo($fname);
//$pathinfo['extension'];
$m=explode(".",$fname);
return end($m);
}
function mainname($fname){
$pathinfo=my_pathinfo($fname);
//print_r($pathinfo);
//exit();
return $pathinfo['filename'];
}
$fp = glob("*.grd");
// grd -> csv
$max_k = count($fp);
foreach($fp as $k=> $v){
$bn = basename($v);
$mn = mainname($v);
/*$csv_outputfile = "{$mn}.csv";
$data = "X,Y,Z\n";
$grd_data = str_replace(" ",",",trim(str_replace("\r","",file_get_contents($v))));
file_put_contents($csv_outputfile,$data.$grd_data);
*/
/*
// 然後建 vrt
$vrt_outputfile = "{$mn}.vrt";
$vrt_template = <<<EOF
<OGRVRTDataSource>
<OGRVRTLayer name="dem">
<SrcDataSource>{$bn}</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>EPSG:3826</LayerSRS>
<GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
</OGRVRTLayer>
</OGRVRTDataSource>
EOF;
file_put_contents($vrt_outputfile,$vrt_template);
*/
// 然後 csv 轉 gpkg
//$CMD = "ogr2ogr -f GPKG {$mn}.gpkg {$mn}.csv -oo AUTODETECT_TYPE=YES -oo X_POSSIBLE_NAMES=field_1 -oo Y_POSSIBLE_NAMES=field_2 -oo Z_POSSIBLE_NAMES=field_3 -a_srs EPSG:3826";
//`{$CMD}`;
// 直接轉 grd 轉 tif
if(!is_file("{$mn}.tif")){
$CMD = "gdal_translate -of GTiff -a_srs EPSG:3826 {$v} {$mn}.tif";
`{$CMD}`;
}
echo sprintf("%d / %d ...\n" , ($k+1),$max_k);
// 用 wrap 重寫一份 #-t_srs EPSG:3826
if(!is_file("{$mn}_fix.tif")){
$CMD = "gdalwarp -r near -dstnodata -32768 {$mn}.tif {$mn}_fix.tif";
`{$CMD}`;
}
}
if(is_file("all.vrt")){
unlink("all.vrt");
}
if(is_file("contour_20m.gpkg")){
unlink("contour_20m.gpkg");
}
if(is_file("8429.db")){
unlink("8429.db");
}
// 合併 成 all.vrt
$CMD = "dir /b *_fix.tif > list.txt";
`{$CMD}`;
$CMD = "gdalbuildvrt -input_file_list list.txt all.vrt";
`{$CMD}`;
// 輸出等高線圖
$CMD = "gdal_contour -i 20 -a elev all.vrt contour_20m.gpkg";
`{$CMD}`;
// 轉出 db
$CMD = "ogr2ogr -f SQLITE 8429_raw.db contour_20m.gpkg -s_srs EPSG:3826 -t_srs EPSG:4326 -nln 8429 -lco GEOMETRY_NAME=GEOMETRY -dialect sqlite -sql \"SELECT ID AS ogc_fid, elev, GEOM AS `GEOMETRY` FROM contour\" ";
`{$CMD}`;
然後 mapfile
MAP
CONFIG PROJ_LIB "/usr/share/proj"
#CONFIG "PROJ_LIB" "proj"
EXTENT 120.457535 23.997763 121.458756 24.448916
FONTSET "/var/www/html/easymap_server/inc/..//fonts/font/fonts.txt"
IMAGETYPE "PNG"
NAME "8429"
MAXSIZE 8192
SIZE 1920 1080
STATUS ON
UNITS METERS
#CONFIG "MS_ERRORFILE" "c:\ms4w_MSSQL\Apache\cgi-bin\easymap_server\log\8429.log"
#debug 5
###GEOTIFFONLY_IMAGECOLOR###
###GEOTIFFONLY_TRANSPARENT###
WEB
METADATA
"wms_srs" "EPSG:3826 EPSG:4326 EPSG:3825 EPSG:3827 EPSG:3828 EPSG:3857 EPSG:900913"
"wms_onlineresource" "https://3wa.tw/easymap_server/wms.php?id=8429&"
"wms_enable_request" "*"
"wms_getcapabilities_version" "1.1.1"
"wms_resx" "1.0"
"wms_resy" "1.0"
"wms_encoding" "UTF-8"
"wms_title" "MOI WMS"
"tile_map_edge_buffer" "10"
"tile_metatile_level" "0"
# WFS Metadata
"wfs_title" "MOI WFS"
"wfs_onlineresource" "https://3wa.tw/easymap_server/wfs.php?id=8429&"
"wfs_srs" "EPSG:4326 EPSG:3826 EPSG:3825 EPSG:3827 EPSG:3828 EPSG:3857 EPSG:900913 CRS:84"
"wfs_enable_request" "*"
END # METADATA
END # WEB
OUTPUTFORMAT
NAME "PNG"
MIMETYPE "image/png"
DRIVER "AGG/PNG"
EXTENSION "png"
IMAGEMODE RGBA
TRANSPARENT TRUE
END # OUTPUTFORMAT
OUTPUTFORMAT
NAME "geojson"
DRIVER "OGR/GEOJSON"
MIMETYPE "application/json; subtype=geojson"
FORMATOPTION "STORAGE=stream"
FORMATOPTION "FORM=SIMPLE"
END
OUTPUTFORMAT
NAME "XML"
DRIVER "OGR/GML"
MIMETYPE "text/xml"
FORMATOPTION "STORAGE=memory"
END # OUTPUTFORMAT
OUTPUTFORMAT
NAME "text/xml"
DRIVER "OGR/GML"
MIMETYPE "text/xml"
FORMATOPTION "STORAGE=memory"
END # OUTPUTFORMAT
PROJECTION
"init=EPSG:4326"
END # PROJECTION
############## CUSTOM START ###############
############## CUSTOM END ###############
############## STYLE START ###############
LEGEND
KEYSIZE 20 10
KEYSPACING 5 5
LABEL
SIZE 10
OFFSET 0 0
SHADOWSIZE 1 1
END # LABEL
STATUS ON
END # LEGEND
QUERYMAP
SIZE -1 -1
STATUS ON
STYLE HILITE
END # QUERYMAP
SCALEBAR
INTERVALS 4
LABEL
SIZE 10
OFFSET 0 0
SHADOWSIZE 1 1
END # LABEL
SIZE 200 3
STATUS OFF
UNITS METERS
END # SCALEBAR
SYMBOL
NAME "circle"
TYPE ellipse
POINTS
1 1
END
END
############## STYLE END ###############
############## 實心圓 START ###############
SYMBOL
NAME "solid_circle"
TYPE ellipse
FILLED true
POINTS 1 1 END
END
############## 實心圓 END ###############
############## 特殊樣式們 START ########################################
############## 空心方塊 START ###############
SYMBOL
NAME "square"
TYPE vector
FILLED false
POINTS
0 0
1 0
1 1
0 1
0 0
END
END
############## 空心方塊 END ###############
############## 實心方塊 START ###############
SYMBOL
NAME "solid_square"
TYPE vector
FILLED true
POINTS
0 0
1 0
1 1
0 1
0 0
END
END
############## 實心方塊 END ###############
############## X字 START ###############
SYMBOL
NAME "kreuz1" # X字形
TYPE VECTOR
POINTS
0 0
0 1
0 -1
0 0
-1 0
1 0
END
END
############## X字 END ###############
############## 十 START ###############
SYMBOL
NAME "kreuz2" # 十
TYPE VECTOR
POINTS
0 0
1 1
-99 -99
0 1
1 0
END
END
############## 十 END ###############
############## dot START ###############
SYMBOL
NAME "dot1" # . 字形
TYPE ellipse
POINTS
1 1
END
FILLED true
ANCHORPOINT 0.5 0.5
END
############## dot END ###############
############## 特殊樣式們 END ########################################
############## POLYLINE START ###############
LAYER
CONNECTION "/var/www/html/easymap_server/inc/..//uploads/8429/8429.db"
CONNECTIONTYPE OGR
DATA "8429"
NAME "lines"
STATUS ON
TYPE LINE
UNITS METERS
# 你要顯示的欄位：改成你真正想顯示的欄位名
LABELITEM "elev"
# 讓線的 label 不要沿線旋轉（固定水平）
# 若你想沿線：把 ANGLE 0 改成 FOLLOW
CLASS
STYLE
ANTIALIAS TRUE
COLOR 0 0 0
WIDTH 0.5
LINECAP round
LINEJOIN round
END
LABEL
TYPE TRUETYPE
FONT "arial"
SIZE 12
COLOR 0 0 0
OUTLINECOLOR 255 255 255
OUTLINEWIDTH 2
# 固定水平顯示（不繞線）
ANGLE 0
# 不要太密：同一條線上重複標註的距離（地圖單位）
REPEATDISTANCE 1200
# 標籤彼此至少距離（地圖單位）
MINDISTANCE 400
# 太短的線不要標
MINFEATURESIZE 80
# 避免切到邊界只剩半個字
PARTIALS FALSE
# 往上偏一點，避免壓在線上
OFFSET 0 0
POSITION AUTO
END
END
# 拉很遠不顯示文字（你可調整門檻）
# 比例尺「大於」此值時(更遠)就不畫 label
LABELMAXSCALEDENOM 200000
METADATA
"wfs_title" "WFS POLYLINE Layer"
"wfs_srs" "EPSG:4326 EPSG:3826 EPSG:3825 EPSG:3827 EPSG:3828 EPSG:3857 EPSG:900913 CRS:84"
"wfs_enable_request" "*"
"wfs_getfeature_formatlist" "GML2,GML3,XML,geojson"
END
END # LAYER
############## POLYLINE END ###############
END # MAP