Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have some tiff files (2 for example) and a shapefile. The real format of the tiff files is geoTiff.

I want to get the points inside the shape boundary. Each point should have 4 parameters: geom(coordinates), red, green, infrared (we have tiff files with such info).

At the end of the operation I want to get a vectorized shape with points.

Could you please give me some advice on how can I do that programmatically. Are there any libs/tools? It is all in the context of geoserver/postgresql/geomajas framework.

Upd: now I have more information. My tiffs will be one of the Landsat combination: http://gif.berkeley.edu/documents/Landsat%20Band%20Information.pdf

The main question now is: to get RGB channels info in exact point, should I parse image points on RGB,or I can get this info from geoTiffMetadata (there is http://download.deegree.org/deegree3/nightly/core/javadoc/org/deegree/coverage/raster/io/imageio/GeoTiffIIOMetadataAdapter.html ). I didn`t find really useful info while debugging the code and calling the methods of GeoTiffIIOMetadataAdapter.

share|improve this question
    
The relation between the image (TIFF) files and the shapefile: shape boundary is situated "inside" the tiff files. The distance between points should be 5-10 meters. –  dementiev May 30 '12 at 21:48
2  
As far as I understand your workflow: 1. Merge rasters 2. Convert your rasters to point features with preserved attributes (x, y, z, R, G, B, I) 3. Extract points within the area of interest (your shape file). I do not know procedure in programs listed by you, but may help in ESRI world. –  Tomek May 31 '12 at 12:25
    
Thank you! I`ll try. –  dementiev May 31 '12 at 15:03
add comment

1 Answer 1

up vote 0 down vote accepted

I tried different approaches(also with postgis db, 3rd party libraries), but stopped at geotools java library. So, I`ll describe how to get points from tiff file that are inside shape polygon boundary. Totally it gets about 15 seconds to do this operation for tiff file of 66 million pixels(200 mb).

public class Run {
    public static void main(String[] args) throws IOException, TransformException, FactoryException, SchemaException {
        long time1 = System.currentTimeMillis();
        Run run = new Run();

        Feature shapeFeature = run.getPolygonFromShape();
        List<MyRGB> myRGBs = run.getOnlyPointsFromShapeWithFragmentation(shapeFeature, run.getTiffMetadata());
        long time2 = System.currentTimeMillis();
        System.out.println("Wasted time = " + (time2 - time1) + " ms");
    }

    private Feature getPolygonFromShape() {
        String filename = "d:\\work\\myShapeFile.shp";
        File file = new File(filename);
        try {
            Map<String, URL> connect = new HashMap<String, URL>();
            connect.put("url", file.toURL());
            DataStore dataStore = DataStoreFinder.getDataStore(connect);
            String[] typeNames = dataStore.getTypeNames();
            String typeName = typeNames[0];
            System.out.println("Reading content " + typeName);
            FeatureSource featureSource = dataStore.getFeatureSource(typeName);
            FeatureCollection collection = featureSource.getFeatures();
            FeatureIterator iterator = collection.features();
            Feature feature = null;
            try {
                while (iterator.hasNext()) {
                    feature = iterator.next();
                }
                return feature;
            } finally {
                iterator.close();
            }

        } catch (Exception e) {
            e.printStackTrace();
        }
        return null;
    }
    public GeoTiffIIOMetadataAdapter getTiffMetadata() throws IOException {
        String filename = "d:\\work\\myTiffFile.tif";
        FileImageInputStream f = new FileImageInputStream(
                new RandomAccessFile(filename, "r"));
        // Look through ImageIO readers
        Iterator iter = ImageIO.getImageReaders(f);
        IIOMetadata imdata = null;
        GeoTiffIIOMetadataAdapter geo_data;
        GeoTiffIIOMetadataAdapter metadata = null;
        while (iter.hasNext() && imdata == null) {
            ImageReader reader = (ImageReader) iter.next();
            reader.setInput(f, true);
            String reader_name = reader.getFormatName().toLowerCase();
            if (reader_name.equalsIgnoreCase("tif")) {
                // Get Image metadata
                imdata = reader.getImageMetadata(0);
                geo_data = new GeoTiffIIOMetadataAdapter(imdata);
                if (geo_data.getGeoKeyDirectoryVersion() == 1) {
                    metadata = geo_data;
                }
            }
        }
        f.close();
        return metadata;
    }

    private List<MyRGB> getOnlyPointsFromShape(Feature shapeFeature, GeoTiffIIOMetadataAdapter tiffMetadata) throws IOException, FactoryException, TransformException {
        BufferedImage bi = ImageIO.read(new File("d:\\work\\myTiffFile.tif"));       
        CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326", true);
        CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:32612");
        MathTransform transformToEpsg4326 = CRS.findMathTransform(targetCRS, sourceCRS);

        double lon = tiffMetadata.getModelTiePoints()[3];// 175784.99999999997 - X
        double lat = tiffMetadata.getModelTiePoints()[4];// 5842215.0   - Y
        GeometryFactory geometryFactory = new GeometryFactory(new PrecisionModel(), 32612);
        Geometry pointUpperLeft = geometryFactory.createPoint(new Coordinate(lon, lat));

        int[] pixel;
        LinkedList<MyRGB> myRGBList = new LinkedList<MyRGB>();
        MyRGB myRGB = null;
        double startLeftUpperTiffX = JTS.transform(pointUpperLeft, transformToEpsg4326).getCoordinate().x;// -115.79200437315563 - X
        double startTiffX = startLeftUpperTiffX;
        double startLeftUpperTiffY = JTS.transform(pointUpperLeft, transformToEpsg4326).getCoordinate().y;// 52.632842905984255   - Y
        double stepX = tiffMetadata.getModelPixelScales()[0];//30.0
        double stepY = tiffMetadata.getModelPixelScales()[1];//30.0
        double radius = 6378137;//earth radius

        MultiPolygon mp = (MultiPolygon) ((SimpleFeatureImpl) shapeFeature).getDefaultGeometry();
        Geometry tmpPoint = null;
        for (int y = 0; y < bi.getHeight(); y++) {
            startLeftUpperTiffX = startTiffX;
            for (int x = 0; x < bi.getWidth(); x++) {
                //create point
                tmpPoint = geometryFactory.createPoint(new Coordinate(startLeftUpperTiffX, startLeftUpperTiffY));
                if (mp.contains(tmpPoint)) {
                    pixel = bi.getRaster().getPixel(x, y, new int[3]);
                    myRGB = new MyRGB(pixel[0], pixel[1], pixel[2]);
                    myRGB.setLatitude(startLeftUpperTiffY);
                    myRGB.setLongitude(startLeftUpperTiffX);
                    myRGBList.add(myRGB);
                }
                double dLon = stepX / (radius * Math.cos(Math.PI * startLeftUpperTiffY / 180));//coordinate, and do plus meters, get new coordinate
                startLeftUpperTiffX += dLon * 180 / Math.PI;
            }
            startLeftUpperTiffY -= (180 / Math.PI) * (stepY / radius);   //have coordinate, and do plus meters, get new coordinate
        }
        System.out.println("myRGBList.size() = " + myRGBList.size());
        System.out.println("startLeftUpperTiffX = " + startLeftUpperTiffX);
        System.out.println("startLeftUpperTiffY = " + startLeftUpperTiffY);
        return myRGBList;
    }

 public class MyRGB {
        int red;
        int green;
        int infrared;
        double latitude;
        double longitude;
//.... constructors, getters and setters
}
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.