Me gustaría crear algunos rásteres de elevación/altura usando python, gdal y numpy. Estoy atascado en numpy (. Y probablemente pitón y gdal)creando el campo de elevación/altura gdal numpy python
En numpy, he estado intentando lo siguiente:
>>> a= numpy.linspace(4,1,4, endpoint=True)
>>> b= numpy.vstack(a)
>>> c=numpy.repeat(b,4,axis=1)
>>> c
array([[ 4., 4., 4., 4.],
[ 3., 3., 3., 3.],
[ 2., 2., 2., 2.],
[ 1., 1., 1., 1.]]) #This is the elevation data I want
de osgeo gdal importación de importación gdalconst *
>>> format = "Terragen"
>>> driver = gdal.GetDriverByName(format)
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE = 1', 'MAXUSERPIXELVALUE= 4'])
>>> raster = numpy.zeros((4,4), dtype=numpy.float32) #this is where I'm messing up
>>> dst_ds.GetRasterBand(1).WriteArray(raster) # gives the null elevation data I asked for in (raster)
0
>>> dst_ds = None
Me imagino que me falta algo simple y espero su consejo.
Gracias,
Chris
(continúa más adelante)
- terragendataset.cpp, v 1,2 *
- Proyecto: YTD (tm) TER conductor
- Propósito: Lector para los documentos Terragen TER
- Autor: Ray Gardener, Daylon Gráficos Ltd. *
- partes de este módulo derivado de los conductores por GDAL
- Frank Warmerdam, ver http://www.gdal.org
Mis disculpas de antemano a Ray Gardener y Frank Warmerdam.
Terragen formato notas:
para la escritura: SCAL = distancia en metros gridpost hv_px = hv_m/SCAL span_px = span_m/SCAL offset = ver TerragenDataset :: write_header() escala = ver TerragenDataset: : write_header() hv = física (hv_px - offset) * 65536.0/escala
le decimos a las personas que llaman que:
Elevations are Int16 when reading,
and Float32 when writing. We need logical
elevations when writing so that we can
encode them with as much precision as possible
when going down to physical 16-bit ints.
Implementing band::SetScale/SetOffset won't work because
it requires callers to know format write details.
So we've added two Create() options that let the
caller tell us the span's logical extent, and with
those two values we can convert to physical pixels.
ds::SetGeoTransform() lets us establish the
size of ground pixels.
ds::SetProjection() lets us establish what
units ground measures are in (also needed
to calc the size of ground pixels).
band::SetUnitType() tells us what units
the given Float32 elevations are in.
Esto me dice que antes de mi WriteArray anterior (somearray) tengo que fijar tanto la transformación geográfica y setProjection y SetUnitType que las cosas funcionen (potencialmente)
Desde el tutorial API GDAL: Python OSR importación importación numpy
dst_ds.SetGeoTransform([ 444720, 30, 0, 3751320, 0, -30 ])
srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS('NAD27')
dst_ds.SetProjection(srs.ExportToWkt())
raster = numpy.zeros((512, 512), dtype=numpy.uint8) #we've seen this before
dst_ds.GetRasterBand(1).WriteArray(raster)
# Once we're done, close properly the dataset
dst_ds = None
Mis disculpas para la creación de un puesto excesivamente largo y una confesión. Si me sale bien, sería bueno tener todas las notas en un solo lugar (la publicación larga).
La Confesión:
he tomado previamente una imagen (JPEG), convertido a una geotiff, y se importa como baldosas en una base de datos PostGIS. Ahora estoy intentando crear un ráster de elevación sobre el cual extender la imagen. Esto probablemente parezca una tontería, pero hay un artista al que quiero honrar, mientras que al no ofenden a los que han trabajado asiduamente para crear y mejorar estas excelentes herramientas.
El artista es belga, por lo que los metros estarían en orden. Ella trabaja en el bajo Manhattan, Nueva York así que, UTM 18. Ahora algo razonable SetGeoTransform. La imagen mencionada anteriormente es w = 3649/h = 2736, tendré que pensar en esto.
Otro intento:
>>> format = "Terragen"
>>> driver = gdal.GetDriverByName(format)
>>> dst_ds = driver.Create('test_3.ter', 4,4,1, gdal.GDT_Float32, ['MINUSERPIXELVALUE=1','MAXUSERPIXELVALUE-4'])
>>> type(dst_ds)
<class 'osgeo.gdal.Dataset'>
>>> import osr
>>> dst_ds.SetGeoTransform([582495, 1, 0.5, 4512717, 0.5, -1]) #x-res 0.5, y_res 0.5 a guess
0
>>> type(dst_ds)
<class 'osgeo.gdal.Dataset'>
>>> srs = osr.SpatialReference()
>>> srs.SetUTM(18,1)
0
>>> srs.SetWellKnownGeogCS('WGS84')
0
>>> dst_ds.SetProjection(srs.ExportToWkt())
0
>>> raster = c.astype(numpy.float32)
>>> dst_ds.GetRasterBand(1).WriteArray(raster)
0
>>> dst_ds = None
>>> from osgeo import gdalinfo
>>> gdalinfo.main(['foo', 'test_3.ter'])
Driver: Terragen/Terragen heightfield
Files: test_3.ter
Size is 4, 4
Coordinate System is:
LOCAL_CS["Terragen world space",
UNIT["m",1]]
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
AREA_OR_POINT=Point
Corner Coordinates:
Upper Left ( 0.0000000, 0.0000000)
Lower Left ( 0.0000000, 4.0000000)
Upper Right ( 4.0000000, 0.0000000)
Lower Right ( 4.0000000, 4.0000000)
Center ( 2.0000000, 2.0000000)
Band 1 Block=4x1 Type=Int16, ColorInterp=Undefined
Unit Type: m
Offset: 2, Scale:7.62939453125e-05
0
Es evidente que cada vez más cerca, pero no está claro si el SetUTM (18,1) fue recogido. ¿Es este un 4x4 en el Hudson River o un Local_CS (sistema de coordenadas)? ¿Qué es un fallo silencioso?
más lectura
// Terragen files aren't really georeferenced, but
// we should get the projection's linear units so
// that we can scale elevations correctly.
// Increase the heightscale until the physical span
// fits within a 16-bit range. The smaller the logical span,
// the more necessary this becomes.
4x4 metros es bastante pequeño lapso de lógica.
Por lo tanto, tal vez esto es tan bueno como se pone. El SetGeoTransform obtiene las unidades correctas, establece la escala y usted tiene su Espacio Mundial Terragen.
Pensamiento final: No programo, pero en cierta medida puedo seguirlo. Dicho esto, yo un poco preguntaba primero si y luego lo miraron los datos como en mi pequeño espacio Mundial YTD (los siguientes gracias a http://www.gis.usu.edu/~chrisg/python/2009/ semana 4):
>>> fn = 'test_3.ter'
>>> ds = gdal.Open(fn, GA_ReadOnly)
>>> band = ds.GetRasterBand(1)
>>> data = band.ReadAsArray(0,0,1,1)
>>> data
array([[26214]], dtype=int16)
>>> data = band.ReadAsArray(0,0,4,4)
>>> data
array([[ 26214, 26214, 26214, 26214],
[ 13107, 13107, 13107, 13107],
[ 0, 0, 0, 0],
[-13107, -13107, -13107, -13107]], dtype=int16)
>>>
Así que esto es gratificante. Imagino la diferencia entre el numpy c utilizado anteriormente y este resultado se aplica a las acciones de aplicar Float16 en este pequeño lapso lógico .
Y de ahí a 'HF2':
>>> src_ds = gdal.Open('test_3.ter')
>>> dst_ds = driver.CreateCopy('test_6.hf2', src_ds, 0)
>>> dst_ds.SetGeoTransform([582495,1,0.5,4512717,0.5,-1])
0
>>> srs = osr.SpatialReference()
>>> srs.SetUTM(18,1)
0
>>> srs.SetWellKnownGeogCS('WGS84')
0
>>> dst_ds.SetProjection(srs.ExportToWkt())
0
>>> dst_ds = None
>>> src_ds = None
>>> gdalinfo.main(['foo','test_6.hf2'])
Driver: HF2/HF2/HFZ heightfield raster
Files: test_6.hf2
test_6.hf2.aux.xml
Size is 4, 4
Coordinate System is:
PROJCS["UTM Zone 18, Northern Hemisphere",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9108"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
VERTICAL_PRECISION=1.000000
Corner Coordinates:
Upper Left ( 0.0000000, 0.0000000) (79d29'19.48"W, 0d 0' 0.01"N)
Lower Left ( 0.0000000, 4.0000000) (79d29'19.48"W, 0d 0' 0.13"N)
Upper Right ( 4.0000000, 0.0000000) (79d29'19.35"W, 0d 0' 0.01"N)
Lower Right ( 4.0000000, 4.0000000) (79d29'19.35"W, 0d 0' 0.13"N)
Center ( 2.0000000, 2.0000000) (79d29'19.41"W, 0d 0' 0.06"N)
Band 1 Block=256x1 Type=Float32, ColorInterp=Undefined
Unit Type: m
0
>>>
casi completamente gratificante, aunque se ve como si estuviera en La Concordia Perú. Así que tengo para descubrir cómo decir, como más al norte, como Nueva York Norte. ¿SetUTM toma un tercer elemento que sugiere "qué tan lejos" al norte o al sur. Parece que ayer me crucé con un gráfico UTM que tenía zonas de etiquetas con letras, algo así como C en el ecuador y quizás T o S en el área de Nueva York.
De hecho, me pareció que la SetGeoTransform fue esencialmente estableciendo la parte superior izquierda norte y este y esto estaba influyendo en el 'cómo extremo norte/sur' parte la de SetUTM. Desactivado a gdal-dev.
Y más tarde aún:
Paddington Bear fue a Perú porque tenía un boleto. Llegué allí porque dije que era donde quería estar. Terragen, trabajando como lo hace, me dio mi espacio en píxeles. Las siguientes llamadas a srs actuaron sobre hf2 (SetUTM), pero el este y el norte se establecieron bajo Terragen, por lo que el UTM 18 se estableció, pero en un cuadro delimitador en el ecuador. Suficientemente bueno.
gdal_translate me llevó a Nueva York. Estoy en Windows así que una línea de comando.y el resultado;
C:\Program Files\GDAL>gdalinfo c:/python27/test_9.hf2
Driver: HF2/HF2/HFZ heightfield raster
Files: c:/python27/test_9.hf2
Size is 4, 4
Coordinate System is:
PROJCS["UTM Zone 18, Northern Hemisphere",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
Origin = (583862.000000000000000,4510893.000000000000000)
Pixel Size = (-1.000000000000000,-1.000000000000000)
Metadata:
VERTICAL_PRECISION=0.010000
Corner Coordinates:
Upper Left ( 583862.000, 4510893.000) (74d 0'24.04"W, 40d44'40.97"N)
Lower Left ( 583862.000, 4510889.000) (74d 0'24.04"W, 40d44'40.84"N)
Upper Right ( 583858.000, 4510893.000) (74d 0'24.21"W, 40d44'40.97"N)
Lower Right ( 583858.000, 4510889.000) (74d 0'24.21"W, 40d44'40.84"N)
Center ( 583860.000, 4510891.000) (74d 0'24.13"W, 40d44'40.91"N)
Band 1 Block=256x1 Type=Float32, ColorInterp=Undefined
Unit Type: m
C:\Program Files\GDAL>
Así que, estamos de vuelta en Nueva York. Probablemente haya mejores formas de abordar todo esto. I tenía que tener un objetivo que aceptara Crear como postulé/improvisé un conjunto de datos de numpy también. Necesito ver otros formatos que permiten crear. Elevación en GeoTiff es una posibilidad (creo)
Mi agradecimiento a todos los comentarios, sugerencias y suaves empujones hacia la lectura adecuada. ¡Hacer mapas en python es divertido!
Chris
¿Cuál es la pregunta? Estás escribiendo ceros en un archivo ... ¿Qué quieres hacer? ¿Escribir ceros no funciona? Si desea escribir la matriz numpy en el archivo, páselo en lugar de la matriz de ceros que está creando. (Puede que necesite pasar 'c.astype (numpy.float32)', si está creando una banda de 32 bits y 'c' es una matriz de 64 bits (que dependerá de la plataforma en la que se encuentre)). –
[El controlador] (http://www.gdal.org/frmt_terragen.html) solo admite 'gdal.GDT_Int16' -no float32 –
@JoeKingston - los ceros estaban funcionando, pero sí quería pasar c como float 32, como Estoy creando mis propios datos de altura. – Chris