OpenDroneMap
https://opendronemap.org/odm/ OpenDroneMap本身是一个集成工具，它通过python 胶水粘合了 Opensfm、ORB_Slam 等三维重建工具。
工具套件: ODM A command line toolkit to process aerial images. Since it’s creation in 2014, it has become the de-facto standard of open source drone image processing. (用于处理航空图像的命令行工具包. 自2014年创建以来，它已成为开源无人机图像处理的事实标准. )
WebODM A user-friendly, extendable application and API for drone image processing. It provides a web interface to ODM with visualization, storage and data analysis functionality. (用于无人机图像处理的用户友好、可扩展的应用程序和API. 它为ODM提供了一个具有可视化、存储和数据分析功能的web界面. )
NodeODM A lightweight REST API to access ODM. It also provides a minimal web interface to access its functions. (用于访问ODM的轻量级REST API. 它还提供了一个最小的web界面来访问其功能. )
支持的格式
Supported Inputs
- JPEGs
- TIFFs (8 and 16bit)
- Ground Control Points
Supported Outputs
- Orthorectified Imagery //正交校正图像 /正射影像图
- Digital Surface Models //数字表面模型
- Digital Terrain Models //数字地形模型
- Textured 3D models //纹理三维模型
- Classified Point Clouds //已分类的点云
WebODM
最简单 易用的 应该是 WebODM了 … 下载页 打包安装版需要收费 ..
自己源码安装 git hub
https://github.com/OpenDroneMap/WebODM/#getting-started
-
Install the following applications (if they are not installed already): Git Docker Docker-compose Python Pip
-
From the Docker Quickstart Terminal or Git Bash (Windows), or from the command line (Mac / Linux), type:
git clone https://github.com/OpenDroneMap/WebODM --config core.autocrlf=input --depth 1
cd WebODM
./webodm.sh start - Open a Web Browser to http://localhost:8000 http://192.168.30.136:8000
Where Are My Files Stored? When using Docker, all processing results are stored in a docker volume and are not available on the host filesystem. There are two specific docker volumes of interest:
Media (called webodm_appmedia): This is where all files related to a project and task are stored.
Postgres DB (called webodm_dbdata): This is what Postgres database uses to store its data.
可以通过参数绑定数据目录
./webodm.sh restart --media-dir /home/user/webodm_data --db-dir /home/user/webodm_db
图片最佳需求 To create a map, press the “Select Images and GCP” button, or drag and drop some images into a project.
You need at least 5 images, but 16-32 is typically the minimum. //至少需5张图像，但通常最少需要16-32张. Images must overlap by 65% or more. Aim for 70-72% //图像必须重叠65%或更多. 目标是70-72% For great 3D, images must overlap by 83% //对于出色的3D，图像必须重叠83% A GCP File is optional, but can increase georeferencing accuracy //GCP文件是可选的，但可以提高地理参考的准确性
installation-and-getting-started
https://docs.opendronemap.org/installation/#installation-and-getting-started
ODM (Docker)
相关文档 https://hackmd.io/@maning/rygH5bh5U
图片超大, 超多的话, 还是用ODM吧…
docker odm 的镜像: docker pull opendronemap/odm
docker odm 的GPU镜像 docker pull opendronemap/odm:gpu
If you have an NVIDIA card, you can test that docker is recognizing the GPU by running:
docker run --rm --gpus all nvidia/cuda:10.0-base nvidia-smi
正射影像图 Creating High Quality Orthophotos
creating-high-quality-orthophotos
Without any parameter tweaks, ODM chooses a good compromise between quality, speed and memory usage. If you want to get higher quality results, you need to tweak some parameters:
--orthophoto-resolution is the resolution of the orthophoto in cm/pixel. Decrease this value for a higher resolution result.
--ignore-gsd is a flag that instructs ODM to skip certain memory and speed optimizations that directly affect the orthophoto. Using this flag will increase runtime and memory usage, but may produce sharper results.
--texturing-data-term should be set to area in forest areas.
--mesh-size should be increased to 300000-600000 and —mesh-octree-depth should be increased to 10-11 in urban areas to recreate better buildings / roofs.
Calibrating the Camera
Camera calibration is a special challenge with commodity cameras. Temperature changes, vibrations, focus, and other factors can affect the derived parameters with substantial effects on resulting data. Automatic or self calibration is possible and desirable with drone flights, but depending on the flight pattern, automatic calibration may not remove all distortion from the resulting products. James and Robson (2014) in their paper Mitigating systematic error in topographic models derived from UAV and ground‐based image networks address how to minimize the distortion from self-calibration.
…
bala bala 一大堆 论述较准相机的问题
Step
看github 的 https://github.com/OpenDroneMap/ODM?tab=readme-ov-file
Quickstart
Docker
The easiest way to run ODM on is via docker. To install docker, see docs.docker.com. Once you have docker installed and working, you can run ODM by placing some images (JPEGs or TIFFs) in a folder named “images” (for example C:\Users\youruser\datasets\project\images or /home/youruser/datasets/project/images) and simply run from a Command Prompt / Terminal:
# Windows
docker run -ti --rm -v c:/Users/youruser/datasets:/datasets opendronemap/odm --project-path /datasets project
# Mac/Linux
docker run -ti --rm -v /home/youruser/datasets:/datasets opendronemap/odm --project-path /datasets projectYou can pass additional parameters by appending them to the command:
docker run -ti --rm -v /datasets:/datasets opendronemap/odm --project-path /datasets project [--additional --parameters --here]For example, to generate a DSM (--dsm) and increase the orthophoto resolution (--orthophoto-resolution 2) :
docker run -ti --rm -v /datasets:/datasets opendronemap/odm --project-path /datasets project --dsm --orthophoto-resolution 2-
拉取 docker 镜像 见上
docker pull opendronemap/odm -
准备数据文件夹
mkdir /home/jk/opendrone/word官方示例的图片数据仓库
图片放到 /home/jk/opendrone/word/images
- 处理图片 + 绑定外部目录
for CPU
cd /home/jk/opendrone/work
docker run -it --rm \
-v "$(pwd)/images:/code/images" \
-v "$(pwd)/output/odm_orthophoto:/code/odm_orthophoto" \
opendronemap/odm- /code/images 待处理的照片
- /code/odm_orthophoto // 对图像文件进行怎么样的处理，根据项目的需要，输入不同的处理指令
-v "$(pwd)/output/odm_meshing:/code/odm_meshing" \ #3D 网面建模
-v "$(pwd)/output/odm_orthophoto:/code/odm_orthophoto" \ #正射影像图
-v "$(pwd)/output/odm_georeferencing:/code/odm_georeferencing" \ # 地理配准后的点云图
-v "$(pwd)/output/odm_texturing:/code/odm_texturing" \ #纹理网面建模其参数
--orthophoto-resolution <float> 正射影像分辨率是正射影像的分辨率，单位为厘米/像素。减小此值可获得更高分辨率的结果。
--orthophoto-png 生成png格式的正射影像
--proj <string> 设置地理坐标系（如WGS84、UTM等）示例: --proj EPSG:4326
--dem-resolution <float> 设置数字高程模型的分辨率（米/像素）
--ignore-gsd 是一个标志，指示 ODM 跳过直接影响正射影像的某些内存和速度优化。使用此标志将增加运行时间和内存使用量，但可能会产生更清晰的结果。
--texturing-data-term 应设置为森林地区的 area。
--mesh-size 应增加到 300000-600000，--mesh-octree-depth 应增加到 10-11，以重建更好的建筑物/屋顶。docker run -it --rm \
-v "$(pwd)/images:/code/images" \
-v "$(pwd)/output/odm_orthophoto0075:/code/odm_orthophoto" \
opendronemap/odm --orthophoto-resolution 0.075
平均每张TIF图4.4MB, 80米高度, 参考大小:
orthophoto-resolution = 0.01, 约3亩面积, 输出文件约 20GB+ orthophoto-resolution = 0.05, 约3亩面积, 输出文件约 2.5GB+ orthophoto-resolution = 0.06, 约3亩面积, 输出文件约 1.6GB+ orthophoto-resolution = 0.065, 约3亩面积, 输出文件约 200MB+ orthophoto-resolution = 0.75, 约3亩面积, 输出文件约 31MB+ orthophoto-resolution = 0.1, 约3亩面积, 输出文件约 20MB+ (由此建议 orthophoto-resolution = 0.70 左右)
for GPU
don’t work?
nvidia-docker run -it --rm \
-v "$(pwd)/images:/code/images" \
-v "$(pwd)/output/odm_orthophoto:/code/odm_orthophoto" \
--gpus all opendronemap/odm:gpu
- 结果
输出数据在:
/home/jk/opendrone/word/output/odm_orthophoto0075
航拍图
nvidia-docker run -it --rm \
-v "$(pwd)/20240712/color_with_geo:/code/images" \
-v "$(pwd)/output/20240712/color_with_geo:/code/odm_orthophoto" \
opendronemap/odm:gpu --orthophoto-resolution 0.075
##############################################
nvidia-docker run -it --rm \
-v "$(pwd)/tif_color_geo:/code/images" \
-v "$(pwd)/output/tif_color_geo:/code/odm_orthophoto" \
opendronemap/odm:gpu --orthophoto-resolution 0.075 --fast-orthophotofor Windows Setup
(https://github.com/OpenDroneMap/ODM#windows-setup)
ODM can be installed natively on Windows. Just download the latest setup from the releases page. After opening the ODM Console you can process datasets by typing:
run C:\Users\youruser\datasets\project [--additional --parameters --here]
# run --help 查看 help
run --project-path E:\content-for-work\odm\for_win\ --fast-orthophoto --feature-type dspsift --min-num-features 20000 --orthophoto-resolution 0.07
run --project-path E:\content-for-work\odm\for_win\ --fast-orthophoto --feature-type dspsift --orthophoto-resolution 0.07# project-path 的文件夹结构是
E:\CONTENT-FOR-WORK\ODM
└─for_win
└─code
├─images #给图片即可
└─opensfm
odm_orthophoto result
odm_orthophoto.tif：- 这通常是最终的正射影像图，它经过了地理校正，可以用于进一步的GIS分析。
- 它可能包含了颜色校正、亮度调整、对比度增强等处理，以得到更清晰和可用的影像。
- 该文件通常是压缩的，以减小文件大小，便于分享和使用。
odm_orthophoto.original.tif：- 这个文件包含了未处理的原始正射影像数据。
- 它没有经过颜色或亮度的调整，保留了原始影像的全部信息。
- 这个文件对于需要原始数据进行分析的用户来说非常有用。
odm_orthophoto_render.tif：- 这个文件是用于渲染的影像，它可能包含了一些额外的处理，如更高级的锐化、色彩平衡或特定效果的增强，以改善视觉效果。
- 它通常用于展示目的，而不是用于精确的分析工作。
所有参数选项 (Options and Flags)
https://docs.opendronemap.org/arguments/
nvidia-docker run -it --rm \
-v "$(pwd)/test:/code/images" \
-v "$(pwd)/output/test2:/code/odm_orthophoto" \
opendronemap/odm:gpu --help usage: run.sh [options] <dataset name>
ODM is a command line toolkit to generate maps, point clouds, 3D models and DEMs from drone, balloon or kite images.
positional arguments:
<dataset name> Name of dataset (i.e subfolder name within project folder). Default: code
optional arguments:
-h, --help show this help message and exit
--project-path <path>
Path to the project folder. Your project folder should contain subfolders for each dataset. Each dataset should have an "images" folder.
--resize-to <integer>
Legacy option (use --feature-quality instead). Resizes images by the largest side for feature extraction purposes only. Set to -1 to disable. This
does not affect the final orthophoto resolution quality and will not resize the original images. Default: 2048
--end-with <string>, -e <string>
End processing at this stage. Can be one of: dataset, split, merge, opensfm, openmvs, odm_filterpoints, odm_meshing, mvs_texturing,
odm_georeferencing, odm_dem, odm_orthophoto, odm_report, odm_postprocess. Default: odm_postprocess
--rerun <string>, -r <string>
Rerun this stage only and stop. Can be one of: dataset, split, merge, opensfm, openmvs, odm_filterpoints, odm_meshing, mvs_texturing,
odm_georeferencing, odm_dem, odm_orthophoto, odm_report, odm_postprocess. Default: None
--rerun-all Permanently delete all previous results and rerun the processing pipeline.
--rerun-from <string>
Rerun processing from this stage. Can be one of: dataset, split, merge, opensfm, openmvs, odm_filterpoints, odm_meshing, mvs_texturing,
odm_georeferencing, odm_dem, odm_orthophoto, odm_report, odm_postprocess. Default: None
--min-num-features <integer>
Minimum number of features to extract per image. More features can be useful for finding more matches between images, potentially allowing the
reconstruction of areas with little overlap or insufficient features. More features also slow down processing. Default: 8000
--feature-type <string>
Choose the algorithm for extracting keypoints and computing descriptors. Can be one of: sift, orb, hahog. Default: sift
--feature-quality <string>
Set feature extraction quality. Higher quality generates better features, but requires more memory and takes longer. Can be one of: ultra, high,
medium, low, lowest. Default: high
--matcher-type <string>
Matcher algorithm, Fast Library for Approximate Nearest Neighbors or Bag of Words. FLANN is slower, but more stable. BOW is faster, but can
sometimes miss valid matches. Can be one of: flann, bow. Default: flann
--matcher-neighbors <integer>
Number of nearest images to pre-match based on GPS exif data. Set to 0 to skip pre-matching. Neighbors works together with Distance parameter, set
both to 0 to not use pre-matching. Default: 8
--matcher-distance <integer>
Distance threshold in meters to find pre-matching images based on GPS exif data. Set both matcher-neighbors and this to 0 to skip pre-matching.
Default: 0
--use-fixed-camera-params
Turn off camera parameter optimization during bundle adjustment. This can be sometimes useful for improving results that exhibit doming/bowling or
when images are taken with a rolling shutter camera. Default: False
--cameras <json> Use the camera parameters computed from another dataset instead of calculating them. Can be specified either as path to a cameras.json file or as a
JSON string representing the contents of a cameras.json file. Default:
--camera-lens <string>
Set a camera projection type. Manually setting a value can help improve geometric undistortion. By default the application tries to determine a
lens type from the images metadata. Can be one of: auto, perspective, brown, fisheye, spherical. Default: auto
--radiometric-calibration <string>
Set the radiometric calibration to perform on images. When processing multispectral and thermal images you should set this option to obtain
reflectance/temperature values (otherwise you will get digital number values). [camera] applies black level, vignetting, row gradient gain/exposure
compensation (if appropriate EXIF tags are found) and computes absolute temperature values. [camera+sun] is experimental, applies all the
corrections of [camera], plus compensates for spectral radiance registered via a downwelling light sensor (DLS) taking in consideration the angle
of the sun. Can be one of: none, camera, camera+sun. Default: none
--max-concurrency <positive integer>
The maximum number of processes to use in various processes. Peak memory requirement is ~1GB per thread and 2 megapixel image resolution. Default:
80
--depthmap-resolution <positive float>
Legacy option (use --pc-quality instead). Controls the density of the point cloud by setting the resolution of the depthmap images. Higher values
take longer to compute but produce denser point clouds. Default: 640
--use-hybrid-bundle-adjustment
Run local bundle adjustment for every image added to the reconstruction and a global adjustment every 100 images. Speeds up reconstruction for very
large datasets. Default: False
--use-3dmesh Use a full 3D mesh to compute the orthophoto instead of a 2.5D mesh. This option is a bit faster and provides similar results in planar areas.
Default: False
--skip-3dmodel Skip generation of a full 3D model. This can save time if you only need 2D results such as orthophotos and DEMs. Default: False
--skip-report Skip generation of PDF report. This can save time if you don't need a report. Default: False
--ignore-gsd Ignore Ground Sampling Distance (GSD). GSD caps the maximum resolution of image outputs and resizes images when necessary, resulting in faster
processing and lower memory usage. Since GSD is an estimate, sometimes ignoring it can result in slightly better image output quality. Default:
False
--mesh-size <positive integer>
The maximum vertex count of the output mesh. Default: 200000
--mesh-octree-depth <integer: 1 <= x <= 14>
Octree depth used in the mesh reconstruction, increase to get more vertices, recommended values are 8-12. Default: 11
--fast-orthophoto Skips dense reconstruction and 3D model generation. It generates an orthophoto directly from the sparse reconstruction. If you just need an
orthophoto and do not need a full 3D model, turn on this option. Default: False
--crop <positive float>
Automatically crop image outputs by creating a smooth buffer around the dataset boundaries, shrinked by N meters. Use 0 to disable cropping.
Default: 3
--boundary <json> GeoJSON polygon limiting the area of the reconstruction. Can be specified either as path to a GeoJSON file or as a JSON string representing the
contents of a GeoJSON file. Default:
--auto-boundary Automatically set a boundary using camera shot locations to limit the area of the reconstruction. This can help remove far away background
artifacts (sky, background landscapes, etc.). See also --boundary. Default: False
--pc-quality <string>
Set point cloud quality. Higher quality generates better, denser point clouds, but requires more memory and takes longer. Each step up in quality
increases processing time roughly by a factor of 4x.Can be one of: ultra, high, medium, low, lowest. Default: medium
--pc-classify Classify the point cloud outputs using a Simple Morphological Filter. You can control the behavior of this option by tweaking the --dem-*
parameters. Default: False
--pc-csv Export the georeferenced point cloud in CSV format. Default: False
--pc-las Export the georeferenced point cloud in LAS format. Default: False
--pc-ept Export the georeferenced point cloud in Entwine Point Tile (EPT) format. Default: False
--pc-filter <positive float>
Filters the point cloud by removing points that deviate more than N standard deviations from the local mean. Set to 0 to disable filtering.
Default: 2.5
--pc-sample <positive float>
Filters the point cloud by keeping only a single point around a radius N (in meters). This can be useful to limit the output resolution of the
point cloud and remove duplicate points. Set to 0 to disable sampling. Default: 0
--pc-tile Reduce the memory usage needed for depthmap fusion by splitting large scenes into tiles. Turn this on if your machine doesn't have much RAM and/or
you've set --pc-quality to high or ultra. Experimental. Default: False
--pc-geometric Improve the accuracy of the point cloud by computing geometrically consistent depthmaps. This increases processing time, but can improve results in
urban scenes. Default: False
--smrf-scalar <positive float>
Simple Morphological Filter elevation scalar parameter. Default: 1.25
--smrf-slope <positive float>
Simple Morphological Filter slope parameter (rise over run). Default: 0.15
--smrf-threshold <positive float>
Simple Morphological Filter elevation threshold parameter (meters). Default: 0.5
--smrf-window <positive float>
Simple Morphological Filter window radius parameter (meters). Default: 18.0
--texturing-data-term <string>
When texturing the 3D mesh, for each triangle, choose to prioritize images with sharp features (gmi) or those that cover the largest area (area).
Default: gmi
--texturing-outlier-removal-type <string>
Type of photometric outlier removal method. Can be one of: none, gauss_clamping, gauss_damping. Default: gauss_clamping
--texturing-skip-global-seam-leveling
Skip normalization of colors across all images. Useful when processing radiometric data. Default: False
--texturing-skip-local-seam-leveling
Skip the blending of colors near seams. Default: False
--texturing-keep-unseen-faces
Keep faces in the mesh that are not seen in any camera. Default: False
--texturing-tone-mapping <string>
Turn on gamma tone mapping or none for no tone mapping. Can be one of none, gamma. Default: none
--gcp <path string> Path to the file containing the ground control points used for georeferencing. The file needs to use the following format: EPSG:<code> or <+proj
definition> geo_x geo_y geo_z im_x im_y image_name [gcp_name] [extra1] [extra2] Default: None
--geo <path string> Path to the image geolocation file containing the camera center coordinates used for georeferencing. Note that omega/phi/kappa are currently not
supported (you can set them to 0). The file needs to use the following format: EPSG:<code> or <+proj definition> image_name geo_x geo_y geo_z
[omega (degrees)] [phi (degrees)] [kappa (degrees)] [horz accuracy (meters)] [vert accuracy (meters)] Default: None
--use-exif Use this tag if you have a GCP File but want to use the EXIF information for georeferencing instead. Default: False
--dtm Use this tag to build a DTM (Digital Terrain Model, ground only) using a simple morphological filter. Check the --dem* and --smrf* parameters for
finer tuning. Default: False
--dsm Use this tag to build a DSM (Digital Surface Model, ground + objects) using a progressive morphological filter. Check the --dem* parameters for
finer tuning. Default: False
--dem-gapfill-steps <positive integer>
Number of steps used to fill areas with gaps. Set to 0 to disable gap filling. Starting with a radius equal to the output resolution, N different
DEMs are generated with progressively bigger radius using the inverse distance weighted (IDW) algorithm and merged together. Remaining gaps are
then merged using nearest neighbor interpolation. Default: 3
--dem-resolution <float>
DSM/DTM resolution in cm / pixel. Note that this value is capped by a ground sampling distance (GSD) estimate. To remove the cap, check --ignore-
gsd also. Default: 5
--dem-decimation <positive integer>
Decimate the points before generating the DEM. 1 is no decimation (full quality). 100 decimates ~99% of the points. Useful for speeding up
generation of DEM results in very large datasets. Default: 1
--dem-euclidean-map Computes an euclidean raster map for each DEM. The map reports the distance from each cell to the nearest NODATA value (before any hole filling
takes place). This can be useful to isolate the areas that have been filled. Default: False
--orthophoto-resolution <float > 0.0>
Orthophoto resolution in cm / pixel. Note that this value is capped by a ground sampling distance (GSD) estimate. To remove the cap, check
--ignore-gsd also. Default: 5
--orthophoto-no-tiled
Set this parameter if you want a striped GeoTIFF. Default: False
--orthophoto-png Set this parameter if you want to generate a PNG rendering of the orthophoto. Default: False
--orthophoto-kmz Set this parameter if you want to generate a Google Earth (KMZ) rendering of the orthophoto. Default: False
--orthophoto-compression <string>
Set the compression to use for orthophotos. Can be one of: JPEG, LZW, PACKBITS, DEFLATE, LZMA, NONE. Default: DEFLATE
--orthophoto-cutline Generates a polygon around the cropping area that cuts the orthophoto around the edges of features. This polygon can be useful for stitching
seamless mosaics with multiple overlapping orthophotos. Default: False
--tiles Generate static tiles for orthophotos and DEMs that are suitable for viewers like Leaflet or OpenLayers. Default: False
--build-overviews Build orthophoto overviews for faster display in programs such as QGIS. Default: False
--cog Create Cloud-Optimized GeoTIFFs instead of normal GeoTIFFs. Default: False
--verbose, -v Print additional messages to the console. Default: False
--copy-to <path> Copy output results to this folder after processing.
--time Generates a benchmark file with runtime info. Default: False
--debug Print debug messages. Default: False
--version Displays version number and exits.
--split <positive integer>
Average number of images per submodel. When splitting a large dataset into smaller submodels, images are grouped into clusters. This value
regulates the number of images that each cluster should have on average. Default: 999999
--split-overlap <positive integer>
Radius of the overlap between submodels. After grouping images into clusters, images that are closer than this radius to a cluster are added to the
cluster. This is done to ensure that neighboring submodels overlap. Default: 150
--split-image-groups <path string>
Path to the image groups file that controls how images should be split into groups. The file needs to use the following format: image_name
group_name Default: None
--sm-cluster <string>
URL to a ClusterODM instance for distributing a split-merge workflow on multiple nodes in parallel. Default: None
--merge <string> Choose what to merge in the merge step in a split dataset. By default all available outputs are merged. Options: all, pointcloud, orthophoto, dem.
Default: all
--force-gps Use images' GPS exif data for reconstruction, even if there are GCPs present.This flag is useful if you have high precision GPS measurements. If'
there are no GCPs, this flag does nothing. Default: False
--gps-accuracy <positive float>
Set a value in meters for the GPS Dilution of Precision (DOP) information for all images. If your images are tagged with high precision GPS
information (RTK), this value will be automatically set accordingly. You can use this option to manually set it in case the reconstruction fails.
Lowering this option can sometimes help control bowling-effects over large areas. Default: 10
--optimize-disk-space
Delete heavy intermediate files to optimize disk space usage. This affects the ability to restart the pipeline from an intermediate stage, but
allows datasets to be processed on machines that don't have sufficient disk space available. Default: False
--pc-rectify Perform ground rectification on the point cloud. This means that wrongly classified ground points will be re-classified and gaps will be filled.
Useful for generating DTMs. Default: False
--primary-band <string>
When processing multispectral datasets, you can specify the name of the primary band that will be used for reconstruction. It's recommended to
choose a band which has sharp details and is in focus. Default: auto
--skip-band-alignment
When processing multispectral datasets, ODM will automatically align the images for each band. If the images have been postprocessed and are
already aligned, use this option. Default: FalseWARNING: No features found in image
在特定场景下，调整一些参数可能有助于改善特征检测效果
- 调整特征检测算法
- 增加特征点数量
- 调整特征匹配参数 匹配比率（Matching Ratio）：降低匹配比率可能允许更多的匹配对，但需要权衡匹配质量。 过滤参数：确保没有过于严格的过滤参数导致特征点被过度筛选。
nvidia-docker run -it --rm \
-v "$(pwd)/20240912_R/550nmR:/code/images" \
-v "$(pwd)/output/20240912_R/550nmR_hahog_16000:/code/odm_orthophoto" \
opendronemap/odm:gpu --orthophoto-resolution 0.075 --feature-type hahog --min-num-features 16000--matcher-type <string> 匹配器算法，用于近似最近邻或词袋的快速库。FLANN 速度较慢，但更稳定。BOW 速度较快，但有时会错过有效匹配。可以是以下之一：flann、bow。默认值：flann
--feature-type <string>
Choose the algorithm for extracting keypoints and computing descriptors. Can be one of: sift, orb, hahog. Default: sift
"""
SIFT 提供高精度特征检测，适用于复杂场景；ORB 速度快，适合实时应用，但精度较低；HAHOG 是基于 SIFT 的变种，计算量更小，适用于处理大规模数据时的快速匹配。
"""
--min-num-features <integer>
每幅图像提取的最小特征数。特征越多，越有助于找到图像之间的更多匹配，从而可能重建重叠较少或特征不足的区域。特征越多也会减慢处理速度。默认值：8000
--matcher-neighbors <integer> 根据 GPS exif 数据进行预匹配的最近图像数量。设置为 0 可跳过预匹配。Neighbors 与 Distance 参数配合使用，将两者均设置为 0 则不使用预匹配。默认值：8
--matcher-distance <integer> 基于 GPS exif 数据查找预匹配图像的距离阈值（以米为单位）。将 matcher-neighbors 和此项都设置为 0 可跳过预匹配。
默认值：0
--orthophoto-resolution <float > 0. 0> 正射影像分辨率（厘米/像素）。请注意，此值受地面采样距离 (GSD) 估计值限制。要取消上限，请同时选中 --ignore-gsd。默认值：5Obj2Tiles 处理工具
将航拍的3D模型转换为 3dtile https://community.cesium.com/t/obj-to-3d-tiles/31567
we have around 180 GB of 3D mesh in OBJ format and would like to visualize it in CesiumJS. Can you suggest the most efficient way of converting this data to 3D Tiles and visualize in CesiumJS? The options that we are aware of are:
using desktop software such as FME (Safe Software) and SURE (nFrames)
Python packages such as Obj2Tiles (GitHub - OpenDroneMap/Obj2Tiles: Converts OBJ files to OGC 3D tiles by performing splitting, decimation and conversion) and objTo3d-tiles (GitHub - PrincessGod/objTo3d-tiles: Convert obj model file to 3d tiles)
Thanks a lot
Obj2Tiles is a fully fledged tool to convert OBJ files to 3D Tiles format. It creates multiple LODs, splits the mesh and repacks the textures.
ODM 开源的, 装了ODM就有这个工具了
文件目录
└─20241029_obj
│ metadata.xml
│
├─BlockABA
│ BlockABA.mtl
│ BlockABA.obj
│ BlockABA_0_0.jpg
│ BlockABA_0_1.jpg
│ BlockABA_0_2.jpg
│ BlockABA_0_3.jpg
│
Basic usage (using defaults)
It runs all the pipeline stages and generates the tileset.json file in the output folder.
Obj2Tiles BlockABA.obj ./output
Decimation
Stop the pipeline at the decimation stage and generate 8 LODs 将大模型分为8个粗细的 LOD
Obj2Tiles --stage Decimation --lods 8 BlockABA.obj ./output
Splitting
Stop the pipeline at the splitting stage and generate 3 divisions per axis 沿每个轴分割为3个块( 3^3 = 3 * 3 * 3 = 9)
Obj2Tiles --stage Splitting --divisions 3 BlockABA.obj ./outputSplitting
Full pipeline
Run all the pipeline stages and generate the tileset.json file in the output folder.
Obj2Tiles --lods 6 --divisions 5 --lat 22.81010706 --lon 113.94045918 --alt 120 BlockABA.obj ./3dTiles_6L_5D
Obj2Tiles --lods 6 --divisions 3 --lat 22.81010706 --lon 113.94045918 --alt 120 BlockABA_c0.3.obj ./3dTiles_6L_3D
Obj2Tiles --lods 6 --divisions 3 --lat 22.81010706 --lon 113.94045918 --alt 120 BlockABA_c005.obj ./3dTiles_6L_3D
http://192.168.20.130/3d/outputTiles/tileset.json
ODM 源码调试
ODM 是如何读取TIF的地理信息?
-
入口代码在:
./run.py,stages\odm_app.py -
看源码: https://github.com/OpenDroneMap/ODM/blob/65fad52506a7c3b07dcac34430a254ccf879e49c/opendm/location.py#L24 读取不到时, 报的一个 warn
-
找到他的测试用例会快一点:
.\tests\test_types.py::test_reconstruction -
图片对象定义在:
opendm\photo.py加载解析便是在parse_exif_values(self, _path_file)中
def parse_exif_values(self, _path_file):
# Disable exifread log
logging.getLogger('exifread').setLevel(logging.CRITICAL)
try:
self.width, self.height = get_image_size.get_image_size(_path_file)
except Exception as e:
raise PhotoCorruptedException(str(e))
tags = {}
xtags = {}
with open(_path_file, 'rb') as f:
# 这里使用的是 exifread 这个库
tags = exifread.process_file(f, details=True, extract_thumbnail=False)
try:
if 'Image Make' in tags:
try:
self.camera_make = tags['Image Make'].values
self.camera_make = self.camera_make.strip()
except UnicodeDecodeError:
log.ODM_WARNING("EXIF Image Make might be corrupted")
self.camera_make = "unknown"
if 'Image Model' in tags:
try:
self.camera_model = tags['Image Model'].values
self.camera_model = self.camera_model.strip()
except UnicodeDecodeError:
log.ODM_WARNING("EXIF Image Model might be corrupted")
self.camera_model = "unknown"
if 'GPS GPSAltitude' in tags:
self.altitude = self.float_value(tags['GPS GPSAltitude'])
if 'GPS GPSAltitudeRef' in tags and self.int_value(tags['GPS GPSAltitudeRef']) is not None and self.int_value(tags['GPS GPSAltitudeRef']) > 0:
self.altitude *= -1
if 'GPS GPSLatitude' in tags and 'GPS GPSLatitudeRef' in tags:
self.latitude = self.dms_to_decimal(tags['GPS GPSLatitude'], tags['GPS GPSLatitudeRef'])
elif 'GPS GPSLatitude' in tags:
log.ODM_WARNING("GPS position for %s might be incorrect, GPSLatitudeRef tag is missing (assuming N)" % self.filename)
self.latitude = self.dms_to_decimal(tags['GPS GPSLatitude'], GPSRefMock('N'))
if 'GPS GPSLongitude' in tags and 'GPS GPSLongitudeRef' in tags:
self.longitude = self.dms_to_decimal(tags['GPS GPSLongitude'], tags['GPS GPSLongitudeRef'])
elif 'GPS GPSLongitude' in tags:
log.ODM_WARNING("GPS position for %s might be incorrect, GPSLongitudeRef tag is missing (assuming E)" % self.filename)
self.longitude = self.dms_to_decimal(tags['GPS GPSLongitude'], GPSRefMock('E'))
if 'Image Orientation' in tags:
self.orientation = self.int_value(tags['Image Orientation'])
except (IndexError, ValueError) as e:
log.ODM_WARNING("Cannot read basic EXIF tags for %s: %s" % (self.filename, str(e)))
try:
if 'Image Tag 0xC61A' in tags:
self.black_level = self.list_values(tags['Image Tag 0xC61A'])
elif 'BlackLevel' in tags:
self.black_level = self.list_values(tags['BlackLevel'])
elif 'Image BlackLevel' in tags:
self.black_level = self.list_values(tags['Image BlackLevel'])
if 'EXIF ExposureTime' in tags:
self.exposure_time = self.float_value(tags['EXIF ExposureTime'])
if 'EXIF FNumber' in tags:
self.fnumber = self.float_value(tags['EXIF FNumber'])
if 'EXIF ISOSpeed' in tags:
self.iso_speed = self.int_value(tags['EXIF ISOSpeed'])
elif 'EXIF PhotographicSensitivity' in tags:
self.iso_speed = self.int_value(tags['EXIF PhotographicSensitivity'])
elif 'EXIF ISOSpeedRatings' in tags:
self.iso_speed = self.int_value(tags['EXIF ISOSpeedRatings'])
if 'Image BitsPerSample' in tags:
self.bits_per_sample = self.int_value(tags['Image BitsPerSample'])
if 'EXIF DateTimeOriginal' in tags:
str_time = tags['EXIF DateTimeOriginal'].values
utc_time = datetime.strptime(str_time, "%Y:%m:%d %H:%M:%S")
subsec = 0
if 'EXIF SubSecTime' in tags:
subsec = self.int_value(tags['EXIF SubSecTime'])
negative = 1.0
if subsec < 0:
negative = -1.0
subsec *= -1.0
subsec = float('0.{}'.format(int(subsec)))
subsec *= negative
ms = subsec * 1e3
utc_time += timedelta(milliseconds = ms)
timezone = pytz.timezone('UTC')
epoch = timezone.localize(datetime.utcfromtimestamp(0))
self.utc_time = (timezone.localize(utc_time) - epoch).total_seconds() * 1000.0
if 'MakerNote SpeedX' in tags and \
'MakerNote SpeedY' in tags and \
'MakerNote SpeedZ' in tags:
self.speed_x = self.float_value(tags['MakerNote SpeedX'])
self.speed_y = self.float_value(tags['MakerNote SpeedY'])
self.speed_z = self.float_value(tags['MakerNote SpeedZ'])
if 'EXIF ExifImageWidth' in tags and \
'EXIF ExifImageLength' in tags:
self.exif_width = self.int_value(tags['EXIF ExifImageWidth'])
self.exif_height = self.int_value(tags['EXIF ExifImageLength'])
except Exception as e:
log.ODM_WARNING("Cannot read extended EXIF tags for %s: %s" % (self.filename, str(e)))
# Warn if GPS coordinates are suspiciously wrong
if self.latitude is not None and self.latitude == 0 and \
self.longitude is not None and self.longitude == 0:
log.ODM_WARNING("%s has GPS position (0,0), possibly corrupted" % self.filename)
# Extract XMP tags
f.seek(0)
xmp = self.get_xmp(f)
## 忽略后面, 且复杂了, 兼容 大疆. Phantom 4 RTK.. 之类设备拍的照片
.....
- so 往这个tif 写 ‘GPS GPSLatitude’ in tags and ‘GPS GPSLatitudeRef’ 这几个标签即可
写入GPS
但, GDAL 可以用于为 GeoTIFF 文件添加地理坐标信息，这些信息通常包含在 GeoTIFF 专有的(GDALMetadata)标签中，并不被 exifread 识别 , 操蛋!
需要使用 piexif 之类的写入规范的 GPS 信息.
def jpg_to_tif_and_copy_GeoTransform(jpg_file_path,ref_tif_file_path, ouput_tif_file_path):
if not(os.path.exists(ref_tif_file_path)):
print(f"{ref_tif_file_path} 参考TIF文件不存在")
return
if not(os.path.exists(jpg_file_path)):
print(f"{jpg_file_path} 不存在")
return
with Image.open(jpg_file_path) as img:
# 将JPG 数据转换为numpy数组
img_array = np.array(img)
# 地理参考的TIFF文件
ref_tif = gdal.Open(ref_tif_file_path)
geotransform = ref_tif.GetGeoTransform()
projection = ref_tif.GetProjection() # 是空的
# srs = osr.SpatialReference()
# srs.ImportFromEPSG(4326)
# projection = srs.ExportToWkt()
# 获取 JPG 图像的尺寸
rows, cols, bands = img_array.shape
# 创建一个新的 TIFF 文件
driver = gdal.GetDriverByName('GTiff')
out_tif = driver.Create(ouput_tif_file_path, cols, rows, bands, gdal.GDT_Byte)
# 复制地理信息到新的 TIFF 文件
out_tif.SetGeoTransform(geotransform)
out_tif.SetProjection(projection)
# 复制 GDALMetadata; ODM没用到, 但多一份也不嫌多
# copyAndWriteEXIF(ref_tif, out_tif)
# 写入 数据到新的 TIFF 文件
for i in range(bands):
out_tif.GetRasterBand(i+1).WriteArray(img_array[:, :, i])
out_tif.FlushCache()
del out_tif
## write gps tags; ODM用的就是这个
copyGPSWriteTags(ref_tif, ouput_tif_file_path)
print(f"复制 {ref_tif_file_path} GPS 到 {ouput_tif_file_path} ")
"""
复制 写入 TAGS 中
pip install piexif
"""
def copyGPSWriteTags(ref_dataset, tif_file_path):
# 创建 EXIF 数据字典
gps_dict = {
# piexif.GPSIFD.GPSLatitudeRef: 'N',
# piexif.GPSIFD.GPSLatitude: [(23, 1), (0, 1), (0, 1)], # 23度
# piexif.GPSIFD.GPSLongitudeRef: 'E',
# piexif.GPSIFD.GPSLongitude: [(113, 1), (0, 1), (0, 1)] # 113度
}
# 参考数据
metadata = ref_dataset.GetMetadata("EXIF")
if isinstance(metadata, dict):
"""
EXIF_GPSAltitude: (140.542)
EXIF_GPSAltitudeRef: 0x00
EXIF_GPSDOP: (1.1)
EXIF_GPSLatitude: (22) (48) (44.2292)
EXIF_GPSLatitudeRef: N
EXIF_GPSLongitude: (113) (56) (14.7096)
EXIF_GPSLongitudeRef: E
EXIF_GPSSpeed: (11.7751)
EXIF_GPSSpeedRef: K
EXIF_GPSVersionID: 0x02 0x02 0x00 0x00
"""
for key, value in metadata.items():
value = value.strip('()')
if key == "EXIF_GPSAltitude":
gps_dict[piexif.GPSIFD.GPSAltitude] = (int(Fraction(value)), 1000)
elif key == "EXIF_GPSAltitudeRef":
gps_dict[piexif.GPSIFD.GPSAltitudeRef] = int(value, 16)
elif key == "EXIF_GPSDOP":
gps_dict[piexif.GPSIFD.GPSDOP] = (int(Fraction(value)), 10)
elif key == "EXIF_GPSLatitude":
lat_values = [float(v) for v in value.replace('(', '').replace(')', '').split()]
gps_dict[piexif.GPSIFD.GPSLatitude] = [
(int(lat_values[0]), 1), # 度
(int(lat_values[1]), 1), # 分
(int(lat_values[2] * 10000), 10000) # 秒
]
elif key == "EXIF_GPSLatitudeRef":
gps_dict[piexif.GPSIFD.GPSLatitudeRef] = value
elif key == "EXIF_GPSLongitude":
lon_values = [float(v) for v in value.replace('(', '').replace(')', '').split()]
gps_dict[piexif.GPSIFD.GPSLongitude] = [
(int(lon_values[0]), 1), # 度
(int(lon_values[1]), 1), # 分
(int(lon_values[2] * 10000), 10000) # 秒
]
elif key == "EXIF_GPSLongitudeRef":
gps_dict[piexif.GPSIFD.GPSLongitudeRef] = value
elif key == "EXIF_GPSSpeed":
gps_dict[piexif.GPSIFD.GPSSpeed] = (int(Fraction(value)), 1000)
elif key == "EXIF_GPSSpeedRef":
gps_dict[piexif.GPSIFD.GPSSpeedRef] = value
elif key == "EXIF_GPSVersionID":
gps_dict[piexif.GPSIFD.GPSVersionID] = bytes([int(x, 16) for x in value.split()])
img = Image.open(tif_file_path)
# 插入 EXIF 数据到图像
exif_bytes = piexif.dump({"GPS": gps_dict})
img.save(tif_file_path, "TIFF", exif=exif_bytes)
"""
1. 读取 jpg_directory 目录的(jpg) RGB图 转换为 tif 输出到 output_directory
2. 读取 tefp_directory 目录的(tif) 扩展geo数据 写入到转换后的 tif
"""
def entry():
jpg_directory = r"E:\content-for-work\2023-0109三普\20240912航拍\_merge\Color"
tefp_directory = r"E:\content-for-work\2023-0109三普\20240912航拍\_merge\850nm"
output_directory = r"E:\content-for-work\2023-0109三普\20240912航拍\_merge\tif_color_geo"
count = 0
if not os.path.exists(output_directory):
os.makedirs(output_directory)
for root, dirs, files in os.walk(jpg_directory):
for jpg_file_name in files:
if jpg_file_name.endswith(".jpg"):# MAX_0001_Color_D.jpg
file_path = os.path.join(root, jpg_file_name)
count += 1
print(f"count={count}, file_path = {file_path}, name={jpg_file_name}")
tefp_file_name = jpg_file_name.replace("Color", "850nm").replace(".jpg", ".tif") # MAX_0001_850nm_D.tif
output_file_name = jpg_file_name.replace(".jpg", ".tif")# MAX_0001_Color_D.tif
# 基于 gdal
jpg_to_tif_and_copy_GeoTransform(
os.path.join(jpg_directory, jpg_file_name)
, os.path.join(tefp_directory, tefp_file_name)
, os.path.join(output_directory, output_file_name) ) ODM 是怎么探测多相机?
- 图片数据对象定义在:
opendm\photo.py, 它有一个band_name用于区分不同的波段 - 加载解析同样是在
parse_exif_values(self, _path_file)中
class ODM_Photo:
"""ODMPhoto - a class for ODMPhotos"""
def __init__(self, path_file):
self.filename = os.path.basename(path_file)
self.mask = None
# Standard tags (virtually all photos have these)
self.width = None
self.height = None
self.camera_make = ''
self.camera_model = ''
self.orientation = 1
# Geo tags
self.latitude = None
self.longitude = None
self.altitude = None
# Multi-band fields
self.band_name = 'RGB'
self.band_index = 0
def parse_exif_values(self, _path_file):
# Disable exifread log
logging.getLogger('exifread').setLevel(logging.CRITICAL)
try:
self.width, self.height = get_image_size.get_image_size(_path_file)
except Exception as e:
raise PhotoCorruptedException(str(e))
tags = {}
xtags = {}
with open(_path_file, 'rb') as f:
......
# Extract XMP tags
f.seek(0)
xmp = self.get_xmp(f)
for xtags in xmp:
try:
# xmp 中的 Camera:BandName 标签
band_name = self.get_xmp_tag(xtags, ['Camera:BandName', '@Camera:BandName', 'FLIR:BandName'])
if band_name is not None:
self.band_name = band_name.replace(" ", "")
so 它的解析是读取 xmp 元数据中 xmp 中的 Camera:BandName 标签值
修改band名称
# 更正目录下tif的 xmp 数据的band
def entry(tif_dir):
if not(os.path.exists(tif_dir)):
print(f"{tif_dir} 目录不存在")
return
count = 0
for root, dirs, files in os.walk(tif_dir):
for name in files:
source_file = os.path.join(root, name)
split_arr = name.split("_")# 1-MAX_0001_Color_D.jpg
sp_0 = split_arr[0] # 1-MAX
sp_1 = split_arr[1] # 0001
wavelength_part = split_arr[2] # 750nm
# flag = split_arr[3] # flag D, R
if(len(split_arr) <= 3):
print(f"{source_file} 已修改 跳过")
continue
# 特定 band 名称; 光谱名需为 NIR{波段}, 可见光名需为 RGB; 不然ODM 给你胡命名
bandname = 'RGB'
if (wavelength_part.lower() == "color"):
pass
else :
bandname = "NIR"+ wavelength_part[:3]
modif_xmp_ext(source_file, bandname)
ouput_file = os.path.join(root, f"{sp_0}_{sp_1}_{wavelength_part}.tif")
os.rename(source_file, ouput_file) # 特定文件名称, ODM 固定以正则的的方式去找文件, 不然报文件找不到.
count+=1
# if(count > 0):
# return
print(f"the end 共处理 {count} 个文件")
# pip install pyexiv2 --proxy=''
from pyexiv2 import Image
def modif_xmp_ext(source_file, Camera_BandName):
img = Image(source_file)
# Extract XMP-metadata from file:
xmp = img.read_xmp()
xmp['Xmp.Camera.BandName'] = Camera_BandName
print(f" 修改 {source_file} BandName 为 {Camera_BandName} ")
# Update XMP-tags within file:
img.modify_xmp(xmp)
# Free resources:
img.close()
passODM 是怎么匹配其他波段文件?
opendm\multispectral.py
关键两个正则表达式:
-
这里使用正则表达式 去匹配
file_regex = re.compile(r"^(.+)[-_]\w+(\.[A-Za-z]{3,4})$") -
使用正则表达式 去匹配其他波段文件
filename_without_band = re.sub(file_regex, "\\1\\2", p.filename)
def compute_band_maps(multi_camera, primary_band):
............
return s2p, p2s
except Exception as e:
# Fallback on filename conventions
log.ODM_WARNING("%s, will use filenames instead" % str(e))
filename_map = {}
s2p = {}
p2s = {}
# 使用正则表达式 去匹配主文件
file_regex = re.compile(r"^(.+)[-_]\w+(\.[A-Za-z]{3,4})$")
for p in primary_band_photos:
filename_without_band = re.sub(file_regex, "\\1\\2", p.filename)
# Quick check
if filename_without_band == p.filename:
raise Exception("Cannot match bands by filename on %s, make sure to name your files [filename]_band[.ext] uniformly." % p.filename)
filename_map[filename_without_band] = p
for band in multi_camera:
photos = band['photos']
for p in photos:
# 使用正则表达式 去匹配其他波段文件
filename_without_band = re.sub(file_regex, "\\1\\2", p.filename)
# Quick check
if filename_without_band == p.filename:
raise Exception("Cannot match bands by filename on %s, make sure to name your files [filename]_band[.ext] uniformly." % p.filename)
s2p[p.filename] = filename_map[filename_without_band]
if band['name'] != band_name:
p2s.setdefault(filename_map[filename_without_band].filename, []).append(p)
return s2p, p2s
TIF 处理
TIF 的元数据
https://blog.csdn.net/ling620/article/details/103731878
- 文件头（File Header）:
- 占用字节数: 8 字节
- 描述: 包含两个字节的字节序（用于标识文件是大端序还是小端序），两个字节的 TIFF 标识符（通常为 “42”），以及 4 字节的 IFD（图像文件目录）偏移量，用于指示第一个 IFD 的位置。
- 图像文件目录（IFD，Image File Directory）:
- 占用字节数: 可变
- 描述: IFD 是一个表，包含了指向实际图像数据和元数据的指针。每个 IFD 条目由 12 字节组成，具体数量取决于图像中包含的标签数量。每个 IFD 开头包含两个字节，用于表示 IFD 中的条目数量，最后 4 字节指向下一个 IFD 的位置（如果有）。
- IFD 条目（IFD Entry）:
- 占用字节数: 每个条目 12 字节
- 描述: 每个 IFD 条目描述图像的一个属性，例如图像宽度、高度、压缩类型、颜色空间等。每个条目包含一个 2 字节的标签标识符，2 字节的数据类型，4 字节的数据项数目，以及 4 字节的值或偏移量。
- 图像数据（Image Data）:
- 占用字节数: 可变
- 描述: 实际的图像数据。它可以是未压缩的像素数据或压缩形式的数据（例如使用 LZW 压缩）。在 GeoTIFF 中，图像数据可以包含一个或多个通道（波段）。
- 地理空间元数据（GeoTIFF Metadata）:
- 占用字节数: 可变
- 描述: GeoTIFF 扩展在 TIFF 中添加了特定的标签，用于存储地理空间信息，例如坐标系、投影、地理范围等。这些元数据通常存储在 IFD 中的自定义标签中。
GeoTiff 主要由 GeoKeyDirctoryTag、GeoDoubleParamsTag、GeoAsciiParamsTag三部分组成。其中GeoKeyDirctoryTag是主要记录信息的部分，而GeoDoubleParamsTag、GeoAsciiParamsTag 可以看作是GeoKeyDirctoryTag的数据字典。
扩展信息
-
EXIF 是最常用的图像元数据格式之一，尤其在数码摄影中。它记录了相机设置（如光圈、快门速度、ISO 等）、拍摄时间、地理位置（如果相机具有 GPS 功能）等信息。 在 TIFF 文件中，EXIF 数据通常作为 IFD（图像文件目录）的一部分，存储在文件头或文件中的特定区域。
-
XMP 是由 Adobe 开发的元数据标准，使用 XML 格式来描述和存储信息。与 EXIF 和 IPTC 不同，XMP 更灵活且扩展性更强，支持更复杂的元数据描述，如版权、描述、标记、分层标签等。
-
IPTC 元数据标准最初由新闻业开发，用于描述和标记图像的版权信息、作者信息、标题、说明、关键词等。它是用于描述图像内容和版权的重要格式。 在 TIFF 文件中，IPTC 数据可以嵌入在文件的任意位置，但通常嵌入在 TIFF 标签中。
元数据存储在 IFD（图像文件目录）中，可以包括 EXIF、IPTC 和 XMP 数据。TIFF 格式的灵活性允许这些数据存储在文件中的任意位置，但通常与图像数据相邻。
Band
Band（波段） 是图像数据的一个组成部分，它代表图像的不同数据层。每个波段包含了同一地点的不同信息，具体内容取决于数据的类型和用途。以下是对波段的详细解释：
- 单波段（Single Band）：
- 包含一个数据层，通常用于灰度图像或单一类型的测量数据，如高程、温度等。
- 在遥感中，单波段图像可能表示某个特定光谱范围内的反射率，例如红外波段。
- 多波段（Multi-Band）：
- 包含多个数据层，每个波段代表不同的信息或光谱范围。
- 在可见光图像中，通常包含三个波段：红（Red）、绿（Green）和蓝（Blue），简称 RGB 图像。
ODM 合成出来的有两个 Band: 第一个Band是灰度信息 0~65535 , 第二Band是透明的 0 | 255;
读写 python
- for gdal
def printTIFMeta(input_tif_file_path):
print(f"======= 文件 {input_tif_file_path} start =======")
# 使用 GDAL 打开TIF文件
dataset = gdal.Open(input_tif_file_path, gdal.GA_ReadOnly)
# 获取文件元数据
file_metadata = dataset.GetMetadata()
# 读取每个波段的数据
num_bands = dataset.RasterCount
print(f"metadata: {file_metadata}")
for i in range(1, num_bands + 1):
band = dataset.GetRasterBand(i)
band_metadata = band.GetMetadata()
print(f" band: {i} = {band_metadata}")
# 获取所有可用域（包括扩展域）
domains = dataset.GetMetadataDomainList()
print("\nAvailable Domains:")
print(domains)
# 遍历并打印每个域的元数据
for domain in domains:
metadata = dataset.GetMetadata(domain)
print(f"\nMetadata for domain '{domain}':")
if isinstance(metadata, dict):
for key, value in metadata.items():
print(f"{key}: {value}")
elif isinstance(metadata, list):
for item in metadata:
print(item)
else:
print(metadata)
# 获取数据的一些元信息，如投影和地理变换
projection = dataset.GetProjection()
geotransform = dataset.GetGeoTransform()
print(f"projection = {projection}")
print(f"geotransform = {geotransform}")
print(f"======= 文件 {input_tif_file_path} end =======")
- for pillow
from PIL import Image
from PIL.TiffTags import TAGS
import os
import re
# ref TAGS
XMP_KEY = 700
def jpg_to_tif_and_copy_xmp(jpg_file_path,xmp_from_tif_file_path, tif_file_path):
if not(os.path.exists(xmp_from_tif_file_path)):
print(f"{xmp_from_tif_file_path} xmp_tif 文件不存在")
return
# JPG文件
with Image.open(jpg_file_path) as img:
# 确保图像模式为RGB，因为某些TIFF格式可能不支持其他模式
if img.mode != 'RGB':
img = img.convert('RGB')
img.save(tif_file_path, 'TIFF')
if not(os.path.exists(tif_file_path)):
print(f"{jpg_file_path} 不能被转换")
return
# 复制 xmp_from_tif_file_path 的 xmp 到转换的TIF
with Image.open(tif_file_path) as tfp:
with Image.open(xmp_from_tif_file_path) as xftfp:
tfp.tag_v2[XMP_KEY] = xftfp.tag_v2[XMP_KEY]# 这里拿到的元数据 是个IDF的目录指针, 基于gdal 可以读取其中的domain信息
tfp.tag_v2.update()
tfp.save(tif_file_path)
print(f"复制 {xmp_from_tif_file_path} XMP 到 {tif_file_path} ")
def entry_test():
jpg_directory = r"E:\content-for-work\2023-0109三普\20240712航拍\color"
tefp_directory = r"E:\content-for-work\2023-0109三普\20240712航拍\850nm"
output_directory = r"E:\content-for-work\2023-0109三普\20240712航拍\color_tif"
jpg_to_tif_and_copy_xmp(jpg_directory+"/MAX_0002_Color_D.jpg"
, tefp_directory+"/MAX_0002_850nm_D.tif"
, output_directory+"/MAX_0002_Color_D.tif")- 写入元数据
XMP信息
def jpg_to_tif_and_copy_ext(jpg_file_path,ref_tif_file_path, ouput_tif_file_path):
if not(os.path.exists(ref_tif_file_path)):
print(f"{ref_tif_file_path} 参考TIF文件不存在")
return
# JPG文件
with Image.open(jpg_file_path) as img:
# 确保图像模式为RGB，因为某些TIFF格式可能不支持其他模式
if img.mode != 'RGB':
img = img.convert('RGB')
img.save(ouput_tif_file_path, 'TIFF')
if not(os.path.exists(ouput_tif_file_path)):
print(f"{ouput_tif_file_path} 不存在")
return
# 复制 的 xmp
copy_xmp_ext(ref_tif_file_path, ouput_tif_file_path)
def copy_xmp_ext(ref_tif_file_path, ouput_tif_file_path):
if not(os.path.exists(ref_tif_file_path)):
print(f"{ref_tif_file_path} 参考TIF文件不存在")
return
if not(os.path.exists(ouput_tif_file_path)):
print(f"{ouput_tif_file_path} 不存在")
return
# 复制 xmp_from_tif_file_path 的 xmp 到转换的TIF
with Image.open(ouput_tif_file_path) as tfp:
with Image.open(ref_tif_file_path) as xftfp:
tfp.tag_v2[XMP_KEY] = xftfp.tag_v2[XMP_KEY]
tfp.tag_v2.update()
tfp.save(ouput_tif_file_path)
print(f"复制 {ref_tif_file_path} XMP 到 {ouput_tif_file_path} ")EXIF 信息
"""
复制 写入 TAGS 中
pip install piexif
"""
def doCopyWriteGPSTags(ref_dataset, tif_file_path):
# 创建 EXIF 数据字典
gps_dict = {
# piexif.GPSIFD.GPSLatitudeRef: 'N',
# piexif.GPSIFD.GPSLatitude: [(23, 1), (0, 1), (0, 1)], # 23度
# piexif.GPSIFD.GPSLongitudeRef: 'E',
# piexif.GPSIFD.GPSLongitude: [(113, 1), (0, 1), (0, 1)] # 113度
}
# 参考数据
metadata = ref_dataset.GetMetadata("EXIF")
if isinstance(metadata, dict):
"""
EXIF_GPSAltitude: (140.542)
EXIF_GPSAltitudeRef: 0x00
EXIF_GPSDOP: (1.1)
EXIF_GPSLatitude: (22) (48) (44.2292)
EXIF_GPSLatitudeRef: N
EXIF_GPSLongitude: (113) (56) (14.7096)
EXIF_GPSLongitudeRef: E
EXIF_GPSSpeed: (11.7751)
EXIF_GPSSpeedRef: K
EXIF_GPSVersionID: 0x02 0x02 0x00 0x00
"""
for key, value in metadata.items():
value = value.strip('()')
if key == "EXIF_GPSAltitude":
gps_dict[piexif.GPSIFD.GPSAltitude] = (int(Fraction(value)), 1000)
elif key == "EXIF_GPSAltitudeRef":
gps_dict[piexif.GPSIFD.GPSAltitudeRef] = int(value, 16)
elif key == "EXIF_GPSDOP":
gps_dict[piexif.GPSIFD.GPSDOP] = (int(Fraction(value)), 10)
elif key == "EXIF_GPSLatitude":
lat_values = [float(v) for v in value.replace('(', '').replace(')', '').split()]
gps_dict[piexif.GPSIFD.GPSLatitude] = [
(int(lat_values[0]), 1), # 度
(int(lat_values[1]), 1), # 分
(int(lat_values[2] * 10000), 10000) # 秒
]
elif key == "EXIF_GPSLatitudeRef":
gps_dict[piexif.GPSIFD.GPSLatitudeRef] = value
elif key == "EXIF_GPSLongitude":
lon_values = [float(v) for v in value.replace('(', '').replace(')', '').split()]
gps_dict[piexif.GPSIFD.GPSLongitude] = [
(int(lon_values[0]), 1), # 度
(int(lon_values[1]), 1), # 分
(int(lon_values[2] * 10000), 10000) # 秒
]
elif key == "EXIF_GPSLongitudeRef":
gps_dict[piexif.GPSIFD.GPSLongitudeRef] = value
elif key == "EXIF_GPSSpeed":
gps_dict[piexif.GPSIFD.GPSSpeed] = (int(Fraction(value)), 1000)
elif key == "EXIF_GPSSpeedRef":
gps_dict[piexif.GPSIFD.GPSSpeedRef] = value
elif key == "EXIF_GPSVersionID":
gps_dict[piexif.GPSIFD.GPSVersionID] = bytes([int(x, 16) for x in value.split()])
img = Image.open(tif_file_path)
# 插入 EXIF 数据到图像
exif_bytes = piexif.dump({"GPS": gps_dict})
img.save(tif_file_path, "TIFF", exif=exif_bytes)多(TIF)图层导出
for ArcMap
打开ArcToolbox 视图: Geoprocessing → ArcToolbox
Arctoolbox → Data Management(数据管理工具) → Raster(栅格) → Raster Dataset (栅格数据集) → Mosaic To New Raster(镶嵌至新栅格)
- Input Rasters 输入多文件
- Output Location 输出目录 Raster Dateset Name with Extension 输出文件名包括扩展名
- Numner of Bands 波段 (跟据原始图,属性可以看到)
- Mosaic Operator (镶嵌运算符) 是最重要的一点，总共有7中运算符
FIRST: 第一个, 根据镶嵌列表中的第一个栅格数据集确定像素值. LAST: 最后一个, 根据叠置的最后一个栅格数据集确定像素值. BLEND: 混合, 方法使用基于距离权重的算法来确定叠置像素的值. MEAN: 平均值, 根据两个叠置的栅格数据集确定平均像素值. MIXIMUM: 最小值, 根据两个叠置的栅格数据集确定较低的像素值. MAXIMUM: 最大值, 方法根据两个叠置的栅格数据集确定较高的像素值. SUM: 总和, 根据叠置栅格数据集确定添加在一起的所有像素的总值.
for GlobalMapper
Global Mapper 是一款地图绘制软件，可将数据（例如：SRTM数据）显示为光栅地图、高程地图、矢量地图，还可以对地图作编辑、转换、打印、记录GPS及利用数据的GIS(地理信息系统)功能
TODO
国内的各地图坐标系转换
GCJ-02 转 WGS84 一定会有偏差的(这本是国标设计的目的)
其一 BD-09, GCJ-02, WGS84 互转 js utils
其二
npm install --save coordtransform
https://www.npmjs.com/package/coordtransform2
//国测局坐标(火星坐标,比如高德地图在用),百度坐标,wgs84坐标(谷歌国外以及绝大部分国外在线地图使用的坐标)
// var coordtransform = require('coordtransform');
import coordtransform from 'coordtransform';
//百度经纬度坐标转国测局坐标
var bd09togcj02 = coordtransform.bd09togcj02(39.915, 116.404);
//国测局坐标转百度经纬度坐标
var gcj02tobd09 = coordtransform.gcj02tobd09(39.915, 116.404);
//wgs84转国测局坐标
var wgs84togcj02 = coordtransform.wgs84togcj02(39.915, 116.404);
//国测局坐标转wgs84坐标
var gcj02towgs84 = coordtransform.gcj02towgs84(39.915, 116.404);
console.log(bd09togcj02);
console.log(gcj02tobd09);
console.log(wgs84togcj02);
console.log(gcj02towgs84);
//result
//bd09togcj02: [ 39.90865673957631, 116.39762729119315 ]
//gcj02tobd09: [ 39.92133699351021, 116.41036949371029 ]
//wgs84togcj02: [ 39.91640428150164, 116.41024449916938 ]
//gcj02towgs84: [ 39.91359571849836, 116.39775550083061 ]WGS84 to GCJ-02
/**
* double latitude = 39.908823;
* double longitude = 116.397470;
* @author yangfh
*/
public static double[] wgs84ToGcj02(double latitude, double longitude) {
if (outOfChina(latitude, longitude)) {
return new double[]{latitude, longitude};
}
double dLat = transformLat(longitude - 105.0, latitude - 35.0);
double dLon = transformLon(longitude - 105.0, latitude - 35.0);
double radLat = latitude / 180.0 * Math.PI;
double magic = Math.sin(radLat);
magic = 1 - EE * magic * magic;
double sqrtMagic = Math.sqrt(magic);
dLat = (dLat * 180.0) / ((A * (1 - EE)) / (magic * sqrtMagic) * Math.PI);
dLon = (dLon * 180.0) / (A / sqrtMagic * Math.cos(radLat) * Math.PI);
double gcjLat = latitude + dLat;
double gcjLon = longitude + dLon;
return new double[]{gcjLat, gcjLon};
}
private static boolean outOfChina(double lat, double lon) {
if (lon < 72.004 || lon > 137.8347) {
return true;
}
if (lat < 0.8293 || lat > 55.8271) {
return true;
}
return false;
}
private static double transformLat(double x, double y) {
double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.sqrt(Math.abs(x));
ret += (20.0 * Math.sin(6.0 * x * Math.PI) + 20.0 * Math.sin(2.0 * x * Math.PI)) * 2.0 / 3.0;
ret += (20.0 * Math.sin(y * Math.PI) + 40.0 * Math.sin(y / 3.0 * Math.PI)) * 2.0 / 3.0;
ret += (160.0 * Math.sin(y / 12.0 * Math.PI) + 320 * Math.sin(y * Math.PI / 30.0)) * 2.0 / 3.0;
return ret;
}
private static double transformLon(double x, double y) {
double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.sqrt(Math.abs(x));
ret += (20.0 * Math.sin(6.0 * x * Math.PI) + 20.0 * Math.sin(2.0 * x * Math.PI)) * 2.0 / 3.0;
ret += (20.0 * Math.sin(x * Math.PI) + 40.0 * Math.sin(x / 3.0 * Math.PI)) * 2.0 / 3.0;
ret += (150.0 * Math.sin(x / 12.0 * Math.PI) + 300.0 * Math.sin(x / 30.0
* Math.PI)) * 2.0 / 3.0;
return ret;
}
public static void main(String[] args) {
double latitude = 39.908823;
double longitude = 116.397470;
double[] convertedCoords = wgs84ToGcj02(latitude, longitude);
double convertedLatitude = convertedCoords[0];
double convertedLongitude = convertedCoords[1];
System.out.println("转换后的坐标：" + convertedLatitude + ", " + convertedLongitude);
}轨迹平滑算法
opencv
opencv 有个 轮廓的多边形逼近 的处理函数 cvApproxPoly() 可以舍近求远用用..
道格拉斯-普克算法(Douglas–Peucker algorithm)
道格拉斯-普克算法(Douglas–Peucker algorithm，亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点，并减少点的数量的一种算法. 该算法的原始类型分别由乌尔斯·拉默(Urs Ramer)于1972年以及大卫·道格拉斯(David Douglas)和托马斯·普克(Thomas Peucker)于1973年提出，并在之后的数十年中由其他学者予以完善
啊, 这个才是专用的算法
public class RamerDouglasPeuckerFilter{
private final double epsilon;
public RamerDouglasPeuckerFilter(double epsilon) {
this.epsilon = epsilon;
}
@Data
@AllArgsConstructor
public static class Coordinate {
private Double lat;
private Double lng;
}
public List<Coordinate> filter(List<Coordinate> points) {
return douglas(points, this.epsilon);
}
private List<Coordinate> douglas(List<Coordinate> points, double epsilon) {
List<Coordinate> result = new ArrayList<>();
// 找到最大阈值点
double maxH = 0;
int index = 0;
int end = points.size();
for (int i = 1; i < end - 1; i++) {
// 计算点到起点和终点组成线段的高
double h = getDistance(points.get(i), points.get(0), points.get(end - 1));
if (h > maxH) {
maxH = h;
index = i;
}
}
// 如果存在最大阈值点，就进行递归遍历出所有最大阈值点
return getCoordinates(points, epsilon, maxH, index, end, result);
}
private List<Coordinate> getCoordinates(List<Coordinate> points, double epsilon,
double maxH, int index, int end, List<Coordinate> result) {
if (maxH > epsilon) {
List<Coordinate> leftPoints = new ArrayList<>();
List<Coordinate> rightPoints = new ArrayList<>();
// 分成两半 继续找比阈值大的
for (int i = 0; i < end; i++) {
if (i < index) {
leftPoints.add(points.get(i));
} else {
rightPoints.add(points.get(i));
}
}
List<Coordinate> leftResult = douglas(leftPoints, epsilon);
List<Coordinate> rightResult = douglas(rightPoints, epsilon);
rightResult.remove(0);
leftResult.addAll(rightResult);
result = leftResult;
} else {
result.add(points.get(0));
result.add(points.get(end - 1));
}
return result;
}
/**
* 计算点到直线的距离
*/
private static double getDistance(Coordinate p, Coordinate s, Coordinate e) {
double AB = distance(s, e);
double CB = distance(p, s);
double CA = distance(p, e);
// 三角形面积
double S = helen(CB, CA, AB);
// 三角形面积 = AB(底) * 高 / 2
// 所以 高 = 2 * 三角形面积 / AB(底)
return 2 * S / AB;
}
/**
* 计算两点之间的距离
*/
private static double distance(Coordinate p1, Coordinate p2) {
double x1 = p1.getLat();
double y1 = p1.getLng();
double x2 = p2.getLat();
double y2 = p2.getLng();
return Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
/**
* 三角形面积
*/
private static double helen(double CB, double CA, double AB) {
double p = (CB + CA + AB) / 2;
return Math.sqrt(p * (p - CB) * (p - CA) * (p - AB));
}
public static void main(String[] args) {
List<Coordinate> points = new ArrayList<>();
points.add(new Coordinate(1D, 2D));
points.add(new Coordinate(2D, 2D));
points.add(new Coordinate(3D, 4D));
points.add(new Coordinate(4D, -5D));
points.add(new Coordinate(5D, 3D));
points.add(new Coordinate(6D, 3D));
points.add(new Coordinate(7D, 5D));
points.add(new Coordinate(8D, 2D));
points.add(new Coordinate(9D, 0D));
points.add(new Coordinate(10D, 9D));
points.add(new Coordinate(11D, 5D));
RamerDouglasPeuckerFilter rdp = new RamerDouglasPeuckerFilter(1);
for (Coordinate p : rdp.filter(points)) {
System.out.print("(" + p.getLat() + "," + p.getLng() + ") ");
}
}
}算法逻辑:
- 在轨迹曲线在曲线首尾两点 A，B 之间连接一条直线 AB，该直线为曲线的弦；
- 遍历曲线上其他所有点，求每个点到直线AB的距离，找到最大距离的点 C，最大距离记为 maxDistance
- 比较该距离 maxDistance 与预先定义的阈值 epsilon 大小，如果 maxDistance < epsilon，则将该直线 AB 作为曲线段的近似，舍去 AB 之间的所有点，曲线段处理完毕；
- 若 maxDistance >= epsilon，则使 C 点将曲线 AB 分为 AC和 CB 两段，并分别对这两段进行(1)~(3)步处理；
- 当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即为原始曲线的路径.
轨迹抽稀之道格拉斯-普克算法 Ramer–Douglas–Peucker algorithm
3dtile osgtb 转b3dm
https://github.com/fanvanzh/3dtiles 3D-Tiles convertion:
Osgb(OpenSceneGraph Binary)to3D-Tiles: convert huge of osgb file to 3D-Tiles.Esri Shapefileto3D-Tiles: convert shapefile to 3D-Tiles.
# from osgb dataset
3dtile.exe -f osgb -i 202409_osgbs -o 202409_3dtile
3dtile.exe -f osgb -i E:\osgb_path -o E:\out_path -c "{\"offset\": 0}"
GDAL 工具库
https://gdal.org/en/stable/download.html#download GDAL（Geospatial Data Abstraction Library）是一个用于读取和写入地理空间数据的开源库。
从GeoJSON创建几何
from osgeo import ogr
geojson = """{"type":"Point","coordinates":[108420.33,753808.59]}"""
point = ogr.CreateGeometryFromJson(geojson)
print("%d,%d" % (point.GetX(), point.GetY()))根据geojson掩膜截图
# 将GeoJSON的几何图形转换为OGR几何对象（递归处理）
def geo_json_geometry_to_ogr(geo_json_geometry):
geometry_type = geo_json_geometry["type"]
coords = geo_json_geometry["coordinates"]
if geometry_type == "Point":
ogr_geom = ogr.Geometry(ogr.wkbPoint)
ogr_geom.SetPoint_2D(0, coords[0], coords[1])
elif geometry_type == "LineString":
ogr_geom = ogr.Geometry(ogr.wkbLineString)
for point in coords:
ogr_geom.AddPoint(point[0], point[1])
elif geometry_type == "Polygon":
ogr_geom = ogr.Geometry(ogr.wkbPolygon)
outer_ring = ogr.Geometry(ogr.wkbLinearRing)
for point in coords[0]:
outer_ring.AddPoint(point[0], point[1])
ogr_geom.AddGeometry(outer_ring)
for inner_ring_coords in coords[1:]:
inner_ring = ogr.Geometry(ogr.wkbLinearRing)
for point in inner_ring_coords:
inner_ring.AddPoint(point[0], point[1])
ogr_geom.AddGeometry(inner_ring)
elif geometry_type == "MultiPoint":
ogr_geom = ogr.Geometry(ogr.wkbMultiPoint)
for point in coords:
sub_point = ogr.Geometry(ogr.wkbPoint)
sub_point.SetPoint_2D(0, point[0], point[1])
ogr_geom.AddGeometry(sub_point)
elif geometry_type == "MultiLineString":
ogr_geom = ogr.Geometry(ogr.wkbMultiLineString)
for line_coords in coords:
sub_line = ogr.Geometry(ogr.wkbLineString)
for point in line_coords:
sub_line.AddPoint(point[0], point[1])
ogr_geom.AddGeometry(sub_line)
elif geometry_type == "MultiPolygon":
ogr_geom = ogr.Geometry(ogr.wkbMultiPolygon)
for poly_coords in coords:
sub_poly = geo_json_geometry_to_ogr({
"type": "Polygon",
"coordinates": poly_coords
})
ogr_geom.AddGeometry(sub_poly)
return ogr_geom
# 将GeoJSON的特征转换为OGR特征对象
def geo_json_feature_to_ogr(geo_json_feature, layer):
feature = ogr.Feature(layer.GetLayerDefn())
# 设置属性
for key, value in geo_json_feature["properties"].items():
feature.SetField(key, value)
# 设置几何图形
geometry = geo_json_feature["geometry"]
ogr_geometry = geo_json_geometry_to_ogr(geometry)
feature.SetGeometry(ogr_geometry)
return feature"""
根据 `geojson_file` geojson文件中的几何掩膜
裁剪 `input_tiff_file` TIF文件中的区域
保存到 `save_path`
分辨率不一致, 一定会导致从采样.
"""
def clip_by_geometry(input_tiff_file, geojson_file, save_path):
if not(os.path.exists(input_tiff_file)):
print(f"{input_tiff_file} 参考TIF文件不存在")
return
tiff_file_name = os.path.basename(input_tiff_file)
with read_geojson_as_ogr_objects(geojson_file) as layer:
numFeatures = layer.GetFeatureCount()
for i in (range(numFeatures)):
feature = layer.GetNextFeature()
id = feature.GetField('Name')
save_name=os.path.join(save_path, f'c_{tiff_file_name}_'+id+'.tif')
select= "Name=" + "'"+ str(id)+ "'"
result = gdal.Warp(
save_name,
input_tiff_file,
format = 'GTiff',
cutlineDSName = geojson_file,# 指定裁剪矢量文件，#shpfile也行
cutlineWhere=select, # 用于裁剪的矢量
cropToCutline = True, # 是否使用cutlineDSName的extent作为输出的界线
resampleAlg = gdal.GRA_NearestNeighbour, # （最近邻插值）来保持原始数值。这种方法对于分类数据（如土地利用类型）尤其有效
- `ear` (最近邻法)
- `bilinear` (双线性插值)
- `cubic` (立方插值)
- `cubic_spline` (立方样条插值)
- `lanczos` (Lanczos 重采样)
- `average` (均值)
# dstNodata = 255 # 输出数据的nodata值
)
# ==========================================
num_bands = result.RasterCount
# 创建一个新的 GeoTIFF 文件
gtiff_driver =gdal.GetDriverByName('GTiff')
# 搞不懂 这里要读一下, 需要的波段才行
for i in range(1, num_bands + 1):
result.GetRasterBand(i).ReadAsArray() #从采样
# TODO 以第一个 band 的数据类型
type_of_data = result.GetRasterBand(1).DataType
# , resampleAlg = gdal.GRA_NearestNeighbour
output_dataset = gtiff_driver.Create(save_name, result.RasterXSize,
result.RasterYSize, num_bands, type_of_data)# gdal.GDT_Byte
# 设置投影和地理变换
output_dataset.SetProjection(result.GetProjection())
output_dataset.SetGeoTransform(result.GetGeoTransform())
# 复制每个波段的数据
for i in range(1, num_bands + 1):
input_band = result.GetRasterBand(i)
output_band = output_dataset.GetRasterBand(i)
output_band.WriteArray(input_band.ReadAsArray())
# data = input_band.ReadAsArray().astype(bytes)
# data = input_band.ReadAsArray().astype(float)
# 写入数据到新的 GeoTIFF 文件
# output_band.WriteArray(data)
# 如果需要，可以设置 nodata 值
output_band.SetNoDataValue(input_band.GetNoDataValue())
# 刷新缓存并关闭数据集
output_dataset.FlushCache()
output_dataset = None
result = None
pass
# __name__ == '__main__'
input_tiff_file = r"E:\content-for-work\2023-0109三普\ndvi分析\ndvi\NDVI_Clip1.tif"
geojson_file = r'E:\content-for-work\2023-0109三普\ndvi分析\all_massif.json'
save_path = r'E:\content-for-work\2023-0109三普\ndvi分析'
clip_by_geometry(input_tiff_file, geojson_file, save_path)
本地切图瓦片(gdal2tiles)
gdal 自带脚本 gdal2tiles.py脚本, 这个脚本可以根据其投影进行坐标切片。
将 大图片TIF 切为 x, y, z 金字塔层级的瓦片, 被本地加载
gdal2tiles.py [-p <profile>] [-r <resampling_method>] [-s <srs>] [-z <zoom_levels>] [-w <webviewer>] [-e] [-a <opacity>] [-v] [-h] <input_raster_file> <output_directory>
- **-p <profile>**：指定切片的配置文件类型，常见的有：
- `mercator`：适用于大多数网络地图应用，采用墨卡托投影，这是最常用的配置。
- `geodetic`：基于大地测量的配置，适用于一些特殊需求的地图展示。
- `raster`：以原始栅格的投影和坐标系统进行切片，保留了原始数据的更多特性。
- **-r <resampling_method>**：指定重采样方法，常用的有：
- `near`：最近邻重采样，速度快，但可能在图像质量上有一定损失，适用于分类数据等。
- `bilinear`：双线性重采样，能提供较好的图像质量，是比较常用的一种方法。
- `cubic`：三次重采样，能提供更高质量的图像，但计算成本相对较高。
- **-s <srs>**：指定输入栅格数据的坐标系统（空间参考系），如果不指定，`gdal2tiles`会尝试自动获取输入数据的坐标系统。
- **-z <zoom_levels>**：指定切片的缩放级别范围。可以用多种方式表示，比如`0-5`表示从缩放级别 0 到 5 进行切片；也可以单独指定一个级别，如`3`表示只在缩放级别 3 进行切片。
- **-w <webviewer>**：指定是否生成一个简单的网页查看器来预览切片后的瓦片。常见的选项有：
- `none`：不生成网页查看器。
- `openlayers`：生成一个基于 OpenLayers 的网页查看器，方便直接在浏览器中查看切片效果。
- `leaflet`：生成一个基于 Leaflet 的网页查看器。
- **-e**：用于排除一些无效的瓦片，比如在边缘部分可能存在一些不完整或不准确的瓦片，通过这个参数可以将其排除。
- **-a <opacity>**：指定切片后瓦片的透明度，取值范围是 0 到 1，0 表示完全透明，1 表示完全不透明。
- **-v**：使命令处于详细模式，会输出更多关于切片过程的信息，有助于了解切片的进展和可能出现的问题。
- **-h**：显示帮助信息，列出所有参数的含义和用法。
- **<input_raster_file>**：要进行切片的输入栅格数据文件的路径，例如`/home/user/data/input.tif`。
- **<output_directory>**：切片后瓦片的输出目录路径，例如`/home/user/output_tiles`。gdal2tiles -p mercator -r bilinear -z 16-20 -w openlayers E:\content-for-work\2023-0109三普\ndvi分析\ndvi\ss_sh_nvdi_rgb_with_nodata.tif E:\content-for-work\2023-0109三普\ndvi分析\ndvi\tiles
# 南山 荔枝园
gdal2tiles -p mercator -r bilinear -z 18-23 -w openlayers E:\content-for-work\2023-0109三普\20240828航拍_荔枝园\sz_ns_lzy.tif E:\content-for-work\2024-12农机智驾APP\map\sz_ns_lzy
# 再来一个高清的
gdal2tiles -p mercator -r bilinear -z 23-23 -w openlayers E:\content-for-work\2023-0109三普\20240828航拍_荔枝园\sz_ns_lzy.tif E:\content-for-work\2024-12农机智驾APP\map\sz_ns_lzy层级参数 -z 16 - 20 , 它还会直接生成本地 openlayer 的 html 加载网页
// 添加点击事件监听器
map.on('click', function(event) {
// 获取点击位置的坐标
var coordinate = event.coordinate;
// 将坐标从地图投影（通常是EPSG:3857）转换为经纬度（EPSG:4326）
var lonLat = ol.proj.toLonLat(coordinate);
// 打印经纬度到控制台
console.log('经度: ' + lonLat[0] + ', 纬度: ' + lonLat[1]);
});TIF压缩(gdal_translate)
gdal_translate -of GTiff -co "COMPRESS=DEFLATE" -co "PREDICTOR=2" -co "ZLEVEL=9" -co "TILED=YES" -co "ALPHA=YES" E:\content-for-work\2023-0109三普\20241010航拍_建图\ss_sh_20241010_1000.tif E:\content-for-work\2023-0109三普\20241010航拍_建图\ss_sh_20241010_1000_d9.tif-co "COMPRESS=DEFLATE"使用DEFLATE压缩算法。-co "PREDICTOR=2"对浮点型数据使用水平差异预测器，对于整数数据也可以尝试PREDICTOR=3。-co "ZLEVEL=9"设置最高级别的压缩。-co "TILED=YES"启用分块存储。-co "ALPHA=YES"添加Alpha通道以支持透明度。
https://docs.geoserver.org/main/en/user/services/wms/outputformats.html
Turf JS工具库
Turf.js 和 JSTS 区别
Turf.js 和 JSTS 都是用于地理空间分析和几何计算的 JavaScript 库，Turf.js 和 JSTS 的主要区别：
-
用途和定位: Turf.js：Turf.js 是一个专注于地理空间分析和地理信息处理的库。它的主要目标是提供简单易用的地理空间函数，可以在前端和后端环境中使用。Turf.js 适用于处理点、线、面等地理空间数据的计算和分析。 JSTS：JSTS（JavaScript Topology Suite）是一个更加全面的几何计算库，它的目标是提供一套完整的几何操作和拓扑关系的函数。JSTS 可以处理复杂的几何操作，例如缓冲区分析、拓扑关系判断、几何分析等。
-
功能和复杂性：
Turf.js：Turf.js 更加轻量级和简单，提供了一些基本的地理空间函数，适合用于处理简单的地理数据。它的设计目标是保持易于使用和快速上手。 JSTS：JSTS 是一个更加复杂和全面的库，提供了丰富的几何操作和拓扑关系的计算功能。它在处理复杂的几何问题和拓扑操作时非常有用，适合用于处理复杂的地理数据和几何分析。
- 依赖和环境： Turf.js：Turf.js 不依赖其他库，可以直接在浏览器或 Node.js 等环境中使用。它的设计目标是使得在前端和后端环境中都能方便地进行地理空间分析。 JSTS：JSTS 依赖于 org.locationtech.jts.io 和 org.locationtech.jts.geom 等 Java 包，通常作为 Java 库被使用。但是，JSTS 也有一个 JavaScript 版本，可以在浏览器或 Node.js 中使用，但需要额外的配置和处理。
综上 Turf.js 更适合用于简单的地理空间分析和处理，而 JSTS 则适合处理复杂的几何操作和拓扑关系计算。根据具体的需求和场景，你可以选择适合的库来处理地理空间数据和几何分析。
安装使用
安装
$ npm install @turf/turf
使用
import * as turf from '@turf/turf'
实例
找 Feature 质心
https://turfjs.org/docs/#centerOfMass //质心 https://turfjs.org/docs/#center //获取要素或要素集合并返回所有要素的绝对中心点。
var path=[113,22.0, 113,22.1]
var t_polygon = turf.polygon([path]);
var centroid = turf.centerOfMass(t_polygon);
const position = centroid.geometry.coordinates;//中心点
几何合并
https://turfjs.org/docs/#union
var poly1 = turf.polygon([[
[-82.574787, 35.594087],
[-82.574787, 35.615581],
[-82.545261, 35.615581],
[-82.545261, 35.594087],
[-82.574787, 35.594087]
]], {"fill": "#0f0"});
var poly2 = turf.polygon([[
[-82.560024, 35.585153],
[-82.560024, 35.602602],
[-82.52964, 35.602602],
[-82.52964, 35.585153],
[-82.560024, 35.585153]
]], {"fill": "#00f"});
var union = turf.union(poly1, poly2);要素合并
https://turfjs.org/docs/#combine Combines a FeatureCollection of Point , LineString , or Polygon features into MultiPoint , MultiLineString , or MultiPolygon features.
var fc = turf.featureCollection([
turf.point([19.026432, 47.49134]),
turf.point([19.074497, 47.509548])
]);
var combined = turf.combine(fc);判断路径是否重叠 (几何重叠)
// [[lng, lat], [lng, lat],[lng, lat]]
function isOverlap(coordinates) {
const lineSegments = [];
// 创建线段
for (let a = 0; a < coordinates.length - 1; a++) {
const lineSegment = turf.lineString([coordinates[a], coordinates[a + 1]]);
lineSegments.push(lineSegment);
}
let lineEnd = turf.lineString([coordinates[coordinates.length-1], coordinates[0]]);
lineSegments.push(lineEnd);
// 判断线段是否相交
for (let i = 0; i < lineSegments.length - 1; i++) {
for (let j = i + 1; j < lineSegments.length; j++) {
const intersectPoints = turf.lineIntersect(lineSegments[i], lineSegments[j]);
if (intersectPoints.features.length > 0) {
// 排除相交于起始点和结束点的情况
const intersectsWithoutEndpoints = intersectPoints.features.filter((point) => {
const pointCoord = point.geometry.coordinates;
return !coordinates.some(coord => coord[0] === pointCoord[0] && coord[1] === pointCoord[1]);
});
if (intersectsWithoutEndpoints.length > 0) {
return true;
}
}
}
}
return false;
}
geometry-api-java
Import/Export operations
From Esri Shape, To Esri Shape From REST Json, To REST Json From GeoJson, To GeoJson From WKT, To WKT From WKB, To WKB