Tengo un archivo ASCII con latitud, longitud y data_val en el siguiente formato.
35-13.643782N, 080-57.190157W, 118.6
...
Tengo un archivo de imagen GeoTiff y puedo verlo fácilmente.
Quiero colocar un "pin" (puede ser un punto / bandera / estrella o lo que sea más fácil) en la imagen en la posición específica de latitud / longitud que se encuentra en el archivo ASCII.
Esto es lo que he logrado hacer hasta ahora:
Mi imagen de origen se ve así:
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
...
Esto es lo que he logrado improvisar en 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
Pero, no puedo entender cómo colocar un "punto rojo" en 35-13.643782N, 080-57.190157W
Tengo que aprender algunos detalles nuevos aquí (nomenclatura sobre SIG).
python
gdal
latitude-longitude
geotiff-tiff
ascii
Brad Walker
fuente
fuente
Respuestas:
Su
gdalinfo
salida muestra que tiene una sola banda GeoTIFF con una tabla de colores (paleta AKA). No puedo ver los valores en esa tabla de colores, por lo que los siguientes comandos convierten la tabla de una sola banda + color en un GeoTIFF RGB de tres bandas. Los comandos también suponen que su archivo ASCII tiene una fila de encabezado y tiene coordenadas en grados decimales, es posible que deba modificar su archivo si no lo tiene.Entradas:
Proceso:
El último comando hace lo siguiente:
gdal_rasterize
Resultado:
fuente
Tienes un buen comienzo.
gdal.CreateCopy
se encargará de la georreferenciación, por lo que no tiene que configurarla por segunda vez utilizando la geotransformación y la proyección.El proceso completo transformará las coordenadas lon / lat en las coordenadas XY de la referencia espacial ráster. Luego, estos cables XY se transformarán en la fila, los índices de col del ráster utilizando la geotransformada inversa. Algún valor de píxel se escribirá en esa posición.
Nota 1:
El comando
gdal.RasterBand.WriteArray(array, xoff, yoff)
opera desde la esquina superior izquierda. En este ejemplo proporciono una matriz 1x1 con el valor 255, de modoxoff
yyoff
son la fila actual, los índices col para la posición lon / lat. Si desea escribir un cuadrado de 3x3, debe ajustarxoff
yyoff
restar 1. También debe asegurarse de que el tipo de datos de la matriz coincida con el del ráster. Como dijiste que quieres un "punto rojo", supongo que hay tres bandas de uint8.fuente