[FreeGeoCZ] pajtn - gdalnumeric.LoadFile

[MICHÁLEK Jan Mgr. <Michalek.Jan AT uhul.cz>]

Ahojte pánové, mám tady, asi triviální prblm, se kterym ale neumim pohnout.

Potřeboval jsem ořezat větąí mnoľství rastrů v tifu podle brambory v ąejpu.

Vygůglil jsem tady todle

http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html

 

Na vzorových datech mi to v osgeo4w běhá uspokojivě. Data, který tim ale chci procedit to neskousne.

 

Problém vzniká tady

------------------------

# Load the source data as a gdalnumeric array

srcArray = gdalnumeric.LoadFile(raster)

 

# Also load as a gdal image to get geotransform (world file) info

srcImage = gdal.Open(raster)

geoTrans = srcImage.GetGeoTransform()

 

# Create an OGR layer from a Field boundary shapefile

field = ogr.Open("%s.shp" % shp)

lyr = field.GetLayer(shp)

poly = lyr.GetNextFeature()

 

# Convert the layer extent to image pixel coordinates

minX, maxX, minY, maxY = lyr.GetExtent()

ulX, ulY = world2Pixel(geoTrans, minX, maxY)

lrX, lrY = world2Pixel(geoTrans, maxX, minY)

 

# Calculate the pixel size of the new image

pxWidth = int(lrX - ulX)

pxHeight = int(lrY - ulY)

 

clip = srcArray[:, ulY:lrY, ulX:lrX]

 

 

Na tom clipu se mě to zrůcá, protoľe pole z původních dat je trojrozměrný a pole z mých dat jen dvou. Kdyľ tam tu první dvojtečku zruąím, tak to doběhne, ale výsledek je podivně pomrąen.

Nepovedlo se mě dohledat, co mi tam dělá tu paseku.

 

Počítám, ľe tady bude někdo, kdo tomu rozumí. Nemohl by mě ten někdo navýst na řeąení.

 

 

Díkec Jelen