GDAL中像素值的纬度和经度
Latitude and Longitude to pixel-value in GDAL
我正在尝试读取一个带存储 0 到 3 之间的值(255 不是 maks)的波段的 GeoTIFF。我的目标是编写一个小程序,它接收 latitude/longitude 和 returns geotiff 中地理坐标处的拟合像素值。
我是从这里下载的:https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/QQHCIK
如果你想看一看...下载 link。 https://drive.google.com/file/d/104De_YQN1V8tSbj2uhIPO0FfowjnInnX/view?usp=sharing
但是,我的代码不起作用。我不能告诉你哪里出了问题,它只是每次都输出错误的像素值。我想要么我的地理坐标到像素的转换是错误的,要么是一个巨大的逻辑错误。
这是我的第一次尝试,它打印出错误的地理坐标值。
此外,它会因某些地理坐标而崩溃,负经度会使它崩溃,因为这会使 pixelY
为负,从而导致光栅方法出现异常。
GdalConfiguration.ConfigureGdal();
var tif = "C:\Users\Lars\Downloads\pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif";
var lat = 51.0;
var lon = 8.3;
using (var image = Gdal.Open(tif, Access.GA_ReadOnly)) {
var bOneBand = image.GetRasterBand(1);
var width = bOneBand.XSize;
var height = bOneBand.YSize;
// Geocoordinate ( lat, lon ) to pixel in the tiff ?
var geoTransform = new double[6];
image.GetGeoTransform(geoTransform);
Gdal.InvGeoTransform(geoTransform, geoTransform);
var bOne = new int[1];
Gdal.ApplyGeoTransform(geoTransform, lat, lon, out var pixelXF, out var pixelYF);
var pixelX = (int)Math.Floor(pixelXF);
var pixelY = (int)Math.Floor(pixelYF);
// Read pixel
bOneBand.ReadRaster(pixelX, pixelY, 1, 1, bOne, 1, 1,0,0);
Console.WriteLine(bOne[0]); // bOne[0] contains wrong data
}
我的第二次尝试如下所示...但也会针对给定坐标输出错误的像素值。它还会因某些地理坐标而崩溃。
GdalConfiguration.ConfigureGdal();
var tif = "C:\Users\Lars\Downloads\pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif";
var lat = 24.7377; // 131.5847
var lon = 131.5847;
using (var image = Gdal.Open(tif, Access.GA_ReadOnly)) {
var bOneBand = image.GetRasterBand(1);
var bOne = new int[1];
// Spatial reference to transform latlng to map coordinates... from python code
var point_srs = new SpatialReference("");
var file_srs = new SpatialReference(image.GetProjection());
point_srs.ImportFromEPSG(4326);
point_srs.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);
var mapCoordinates = new double[2]{ lat, lon};
var transform = new CoordinateTransformation(point_srs, file_srs);
transform.TransformPoint(mapCoordinates);
// Map coordinates to pixel coordinates ?
var geoTransform = new double[6];
image.GetGeoTransform(geoTransform);
Gdal.InvGeoTransform(geoTransform, geoTransform);
Gdal.ApplyGeoTransform(geoTransform, mapCoordinates[0], mapCoordinates[1], out var pixelXF, out var pixelYF);
var pixelX = (int)pixelXF;
var pixelY = (int)pixelYF;
bOneBand.ReadRaster(pixelX, pixelY, 1, 1, bOne, 1, 1, 0,0);
Console.WriteLine(bOne[0]); // bOne[0] contains wrong value
}
我的代码有什么问题,为什么会输出错误的像素值?
任何帮助表示赞赏!
也许您在应该 lon/lat 的地方使用 lat/lon?如果您打印了一些值,例如 mapCoordinates
.
,那将会很有帮助
以下是使用 R 执行此操作的方法,以供比较:
library(terra)
r = rast("pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif")
xy = cbind(c(8.3, 131.5847), c(51.0, 24.7377))
extract(r, xy)
# pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1
#1 9
#2 NA
和 zero-based row/col 个数字
rowColFromCell(r, cellFromXY(r, xy)) - 1
# [,1] [,2]
#[1,] 4364 22596
#[2,] 7515 37390
并且您可以使用 describe
(a.k.a.GDALinfo) 获取一些元数据。例如
d = describe("pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif")
d[50]
# [1] "Band 1 Block=43200x1 Type=Byte, ColorInterp=Gray"
我正在尝试读取一个带存储 0 到 3 之间的值(255 不是 maks)的波段的 GeoTIFF。我的目标是编写一个小程序,它接收 latitude/longitude 和 returns geotiff 中地理坐标处的拟合像素值。
我是从这里下载的:https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/QQHCIK 如果你想看一看...下载 link。 https://drive.google.com/file/d/104De_YQN1V8tSbj2uhIPO0FfowjnInnX/view?usp=sharing
但是,我的代码不起作用。我不能告诉你哪里出了问题,它只是每次都输出错误的像素值。我想要么我的地理坐标到像素的转换是错误的,要么是一个巨大的逻辑错误。
这是我的第一次尝试,它打印出错误的地理坐标值。
此外,它会因某些地理坐标而崩溃,负经度会使它崩溃,因为这会使 pixelY
为负,从而导致光栅方法出现异常。
GdalConfiguration.ConfigureGdal();
var tif = "C:\Users\Lars\Downloads\pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif";
var lat = 51.0;
var lon = 8.3;
using (var image = Gdal.Open(tif, Access.GA_ReadOnly)) {
var bOneBand = image.GetRasterBand(1);
var width = bOneBand.XSize;
var height = bOneBand.YSize;
// Geocoordinate ( lat, lon ) to pixel in the tiff ?
var geoTransform = new double[6];
image.GetGeoTransform(geoTransform);
Gdal.InvGeoTransform(geoTransform, geoTransform);
var bOne = new int[1];
Gdal.ApplyGeoTransform(geoTransform, lat, lon, out var pixelXF, out var pixelYF);
var pixelX = (int)Math.Floor(pixelXF);
var pixelY = (int)Math.Floor(pixelYF);
// Read pixel
bOneBand.ReadRaster(pixelX, pixelY, 1, 1, bOne, 1, 1,0,0);
Console.WriteLine(bOne[0]); // bOne[0] contains wrong data
}
我的第二次尝试如下所示...但也会针对给定坐标输出错误的像素值。它还会因某些地理坐标而崩溃。
GdalConfiguration.ConfigureGdal();
var tif = "C:\Users\Lars\Downloads\pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif";
var lat = 24.7377; // 131.5847
var lon = 131.5847;
using (var image = Gdal.Open(tif, Access.GA_ReadOnly)) {
var bOneBand = image.GetRasterBand(1);
var bOne = new int[1];
// Spatial reference to transform latlng to map coordinates... from python code
var point_srs = new SpatialReference("");
var file_srs = new SpatialReference(image.GetProjection());
point_srs.ImportFromEPSG(4326);
point_srs.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);
var mapCoordinates = new double[2]{ lat, lon};
var transform = new CoordinateTransformation(point_srs, file_srs);
transform.TransformPoint(mapCoordinates);
// Map coordinates to pixel coordinates ?
var geoTransform = new double[6];
image.GetGeoTransform(geoTransform);
Gdal.InvGeoTransform(geoTransform, geoTransform);
Gdal.ApplyGeoTransform(geoTransform, mapCoordinates[0], mapCoordinates[1], out var pixelXF, out var pixelYF);
var pixelX = (int)pixelXF;
var pixelY = (int)pixelYF;
bOneBand.ReadRaster(pixelX, pixelY, 1, 1, bOne, 1, 1, 0,0);
Console.WriteLine(bOne[0]); // bOne[0] contains wrong value
}
我的代码有什么问题,为什么会输出错误的像素值? 任何帮助表示赞赏!
也许您在应该 lon/lat 的地方使用 lat/lon?如果您打印了一些值,例如 mapCoordinates
.
以下是使用 R 执行此操作的方法,以供比较:
library(terra)
r = rast("pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif")
xy = cbind(c(8.3, 131.5847), c(51.0, 24.7377))
extract(r, xy)
# pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1
#1 9
#2 NA
和 zero-based row/col 个数字
rowColFromCell(r, cellFromXY(r, xy)) - 1
# [,1] [,2]
#[1,] 4364 22596
#[2,] 7515 37390
并且您可以使用 describe
(a.k.a.GDALinfo) 获取一些元数据。例如
d = describe("pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif")
d[50]
# [1] "Band 1 Block=43200x1 Type=Byte, ColorInterp=Gray"