Eu tenho um arquivo ASCII com latitude, longitude e data_val no seguinte formato.
35-13.643782N, 080-57.190157W, 118.6
...
Eu tenho um arquivo de imagem GeoTiff e posso visualizá-lo facilmente.
Quero colocar um "alfinete" (pode ser um ponto / bandeira / estrela ou o que for mais fácil) na imagem na posição específica de latitude / longitude encontrada no arquivo ASCII.
Aqui está o que eu consegui fazer até agora:
Minha imagem de origem é assim:
Driver: GTiff/GeoTIFF
Files: /tmp/Charlotte SEC 100.tif
Size is 16867, 12358
Coordinate System is:
PROJCS["Lambert Conformal Conic",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.2572221010042,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4269"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",38.66666666666666],
PARAMETER["standard_parallel_2",33.33333333333334],
PARAMETER["latitude_of_origin",34.11666666666667],
PARAMETER["central_meridian",-78.75],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-365041.822331817995291,240536.419747152860509)
Pixel Size = (42.334586069440391,-42.334898968590878)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_DATETIME=2016:06:24 12:46:45
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
TIFFTAG_XRESOLUTION=300
TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -365041.822, 240536.420) ( 82d48'55.43"W, 36d13' 4.92"N)
Lower Left ( -365041.822, -282638.262) ( 82d35'10.11"W, 31d30'17.00"N)
Upper Right ( 349015.641, 240536.420) ( 74d51'46.40"W, 36d13'26.16"N)
Lower Right ( 349015.641, -282638.262) ( 75d 4'55.60"W, 31d30'36.99"N)
Center ( -8013.091, -21050.921) ( 78d50'12.11"W, 33d55'36.35"N)
Band 1 Block=16867x1 Type=Byte, ColorInterp=Palette
Color Table (RGB with 256 entries)
0: 255,255,255,255
...
Aqui está o que eu consegui juntar em Python:
from osgeo import gdal, osr
src_filename = '/tmp/Charlotte SEC 100.tif'
dst_filename = '/tmp/foo.tiff'
# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)
# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)
# Specify raster location through geotransform array
# (upperleftx, scalex, skewx, upperlefty, skewy, scaley)
# Scale = size of one pixel in units of raster projection
# this example below assumes 100x100
gt = [-365041.822, 100, 0, 240536.420, 0, -100]
# Set location
dst_ds.SetGeoTransform(gt)
# Get raster projection
epsg = 4269 # http://spatialreference.org/ref/sr-org/lambert_conformal_conic_2sp/
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()
# Set projection
dst_ds.SetProjection(dest_wkt)
# Close files
dst_ds = None
src_ds = None
Mas não consigo descobrir como colocar um "ponto vermelho" em 35-13.643782N, 080-57.190157W
Estou tendo que aprender alguns novos detalhes aqui (nomenclatura sobre SIG).
python
gdal
latitude-longitude
geotiff-tiff
ascii
Brad Walker
fonte
fonte
Respostas:
Sua
gdalinfo
saída mostra que você tem uma única banda GeoTIFF com uma tabela de cores (paleta AKA). Como não vejo os valores nessa tabela de cores, os comandos abaixo convertem a tabela de banda única + cor em um GeoTIFF RGB de três bandas. Os comandos também assumem que seu arquivo ASCII possui uma linha de cabeçalho e tem coordenadas em graus decimais; talvez seja necessário modificar o arquivo, se não houver.Entradas:
Processo:
O último comando faz o seguinte:
gdal_rasterize
Resultado:
fonte
Você começou bem.
gdal.CreateCopy
cuidará da georreferenciação, para que você não precise defini-la uma segunda vez usando a geotransformação e a projeção.O processo completo transformará as cordas lon / lat nas coordenadas XY da referência espacial raster. Em seguida, essas cordas XY serão transformadas na linha, índices de colunas da varredura usando a geotransformação inversa. Algum valor de pixel será gravado nessa posição.
Nota 1:
O comando
gdal.RasterBand.WriteArray(array, xoff, yoff)
opera no canto superior esquerdo. Neste exemplo, forneço uma matriz 1x1 com o valor 255, assimxoff
eyoff
são a linha real, índices de col para a posição lon / lat. Se você deseja escrever um quadrado 3x3, precisa ajustarxoff
eyoff
subtraindo 1. Você também deve garantir que o tipo de dados da matriz corresponda ao da varredura. Como você disse que quer um "ponto vermelho", estou assumindo que existem três faixas de uint8.fonte