使用边界盒坐标的GeoTiff文件的地理发行，该文件从JPEG转换

Question

我目前正在使用仅作为JPEG使用的航空影像。我使用带有坐标系EPSG的边界箱获得数据：28992。因此，我认为我可以使用该边界箱的坐标作为接地控制点，将其与GDAL一起进行地进行。但是，该代码有效，但是当我使用rasterio打开它们时，.TIF文件仍未被授予。这是代码：

HRtif=gdal.Translate('orthoHR.tif', 'orthoHR.jpeg', format='GTiff' )
IRtif=gdal.Translate('ortho25IR.tif', 'ortho25IR.jpeg', format='GTiff' )


#georeferencing using bbox coordinates
HRref= gdal.Open('orthoHR.tif', gdal.GA_Update)
IRref= gdal.Open('ortho25IR.tif', gdal.GA_Update)

#set coordinate system
sr = osr.SpatialReference()
sr.ImportFromEPSG(28992)

#GCPS bbox, p1=x1,y1 p2= x1,y2 p3= x2,y1 p4= x2,y2
x1=1633134.884
y1=3669075.728
x2=1652322.639
y2=3684544.511

gcp1= (x1, y1)
gcp2= (x1, y2) 
gcp3= (x2, y1) 
gcp4= (x2, y2)

#which pixels correspond
HRarray=HRref.ReadAsArray()
IRarray=IRref.ReadAsArray()
px1=0
py1=0
px2= HRarray.shape[1]
py2= HRarray.shape[2]

pp1= (px1, py1)
pp2= (px1, py2)
pp3= (px2, py1)
pp4= (px2, px2)

#make GCPS
gcps= [gdal.GCP(gcp1[0],gcp1[1], 0, pp1[0],pp1[1]),
       gdal.GCP(gcp2[0],gcp2[1], 0, pp2[0],pp2[1]),
       gdal.GCP(gcp3[0],gcp3[1], 0, pp3[0],pp3[1]),
       gdal.GCP(gcp4[0],gcp4[1], 0, pp4[0],pp4[1])]

#apply GCPS
HRref.SetGCPs(gcps, sr.ExportToWkt())
IRref.SetGCPs(gcps, sr.ExportToWkt())

HRref=None
IRref= None


HRras = rasterio.open('orthoHR.tif', driver='GTiff', crs='EPSG:28992')
IRras = rasterio.open('ortho25IR.tif', driver='GTiff', crs='EPSG:28992')



print(HRras.crs)
print(IRras.crs)

这是输出：

None
None
/usr/local/lib/python3.7/dist-packages/rasterio/__init__.py:220: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix be returned.
  s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)```

Answer 1

设置原始文件的CRS：

与 rasterio 打开后，不要使用 crs 尝试在设置GCP之前

#apply GCPS
HRref.SetSpatialRef(sr)
IRref.SetSpatialRef(sr)
HRref.SetGCPs(gcps, sr)
IRref.SetGCPs(gcps, sr)

Answer 2

退后一步并从您从哪里获取数据可能会有所帮助。我不确定这是否是公共服务，但例如：

例如，例如WMS红外数据集的终点：

您可以使用该URL在命令行上运行 gdalinfo ，以查看其中的可用层/数据集。不幸的是，使用python的 gdal.info 仅显示“频段”，这无济于事。但是打开主乌尔尔还允许检索类似的子dataset，例如：

from osgeo import gdal

wms_url = "https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?request=GetCapabilities&service=wms"

ds = gdal.OpenEx(wms_url)
sds = ds.GetSubDatasets()
ds = None

现在 sds 包含子数据图：

[('WMS:https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=2016_ortho25IR&SRS=EPSG:28992&BBOX=-2000.0,290000.0,294000.0,630000.0&TRANSPARENT=FALSE',
  'Luchtfoto 2016 Ortho 25cm Infrarood'),
 ('WMS:https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=2017_ortho25IR&SRS=EPSG:28992&BBOX=-2000.0,290000.0,294000.0,630000.0&TRANSPARENT=FALSE',
  'Luchtfoto 2017 Ortho 25cm Infrarood'),
 ('WMS:https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=2018_ortho25IR&SRS=EPSG:28992&BBOX=-2000.0,290000.0,294000.0,630000.0&TRANSPARENT=FALSE',
  'Luchtfoto 2018 Ortho 25cm Infrarood'),
 ('WMS:https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=2019_ortho25IR&SRS=EPSG:28992&BBOX=-2000.0,290000.0,294000.0,630000.0&TRANSPARENT=FALSE',
  'Luchtfoto 2019 Ortho 25cm Infrarood'),
 ('WMS:https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=2020_ortho25IR&SRS=EPSG:28992&BBOX=-2000.0,290000.0,294000.0,630000.0&TRANSPARENT=FALSE',
  'Luchtfoto 2020 Ortho 25cm Infrarood'),
 ('WMS:https://service.pdok.nl/hwh/luchtfotocir/wms/v1_0?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=Actueel_ortho25IR&SRS=EPSG:28992&BBOX=-2000.0,290000.0,294000.0,630000.0&TRANSPARENT=FALSE',
  'Luchtfoto Actueel Ortho 25cm Infrarood')]

如果您探索主URL（ wms_url ），例如在浏览器中，您'将看到其他预测可用，GDAL只是选择/显示第一个预测。在这种情况下，这似乎也最合理。

在这种情况下， sds 中的最后一个条目包含最新的图像（“ actueel”）。使用该下载子集（以粗略的10m分辨率）可以使用：

wms_layer, layer_name = sds[-1]
output_tif = rf"C:\Temp\{layer_name.replace('' ', '_').lower}_10m_subset.tif"

ds = gdal.Translate(
    output_tif, 
    wms_layer, 
    projWin=[109000, 482000, 115000, 476000], 
    xRes=10,
    yRes=10
)
ds = None

将下载数据，并将其存储在已经具有georepercted的Geotiff中。在QGI中可视化它，从QuickMapservices插件上叠加在“ Google Hybrid”层上，看起来像：

（Google Maps肯定到处都不是完美的，因此只需将其用作快速检查）即可。

因此，如果您以类似的方式获得了数据，则可能已经正确地在EPSG上进行了地理位置：28992 GRID。在这种情况下，您应该能够简单地分配正确的元数据。这比使用GCP的GCP优先，这通常仅适用于尚未正确投影的数据。后者可能有可能引入不准确性或触发不必要的重新采样等。

您可以使用gdal.translate将您必须的jpeg转换为geotiff，分配（不转换！）可以使用类似的东西来完成投影/边界盒

ds = gdal.Translate(out_tif, in_jpeg, outputBounds=bbox, outputSRS="EPSG:28992")
ds = None

。方式：

ds = gdal.OpenEx(in_jpeg, gdal.GA_Update)

xsize = ds.RasterXSize
ysize = ds.RasterYSize

ds.SetProjection(srs.ExportToWkt())

# convert boundingbox to geotransform
# gt = [ulx, xres, _, uly, _, yres]

ds.SetGeoTransform(gt)
ds = None