Source code

001/*
002 * Copyright (c) 2015 The Regents of the University of California.
003 * All rights reserved.
004 *
005 * '$Author: crawl $'
006 * '$Date: 2018-10-16 03:23:34 +0000 (Tue, 16 Oct 2018) $' 
007 * '$Revision: 34712 $'
008 * 
009 * Permission is hereby granted, without written agreement and without
010 * license or royalty fees, to use, copy, modify, and distribute this
011 * software and its documentation for any purpose, provided that the above
012 * copyright notice and the following two paragraphs appear in all copies
013 * of this software.
014 *
015 * IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
016 * FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
017 * ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
018 * THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
019 * SUCH DAMAGE.
020 *
021 * THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
022 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
023 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
024 * PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
025 * CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
026 * ENHANCEMENTS, OR MODIFICATIONS.
027 *
028 */
029package org.kepler.gis.util;
030
031import java.io.BufferedReader;
032import java.io.File;
033import java.io.IOException;
034import java.io.InputStream;
035import java.io.InputStreamReader;
036import java.util.Collections;
037import java.util.HashSet;
038import java.util.Set;
039import java.util.regex.Matcher;
040import java.util.regex.Pattern;
041import java.util.List;
042import java.util.ArrayList;
043
044import org.apache.commons.io.FilenameUtils;
045import org.geotools.coverage.grid.io.AbstractGridCoverage2DReader;
046import org.geotools.coverage.grid.io.AbstractGridFormat;
047import org.geotools.coverage.grid.io.GridFormatFinder;
048import org.geotools.coverage.grid.io.UnknownFormat;
049import org.geotools.factory.Hints;
050import org.geotools.geometry.GeneralEnvelope;
051import org.geotools.geometry.jts.ReferencedEnvelope;
052import org.geotools.referencing.CRS;
053import org.kepler.gis.data.RasterToken;
054import org.opengis.coverage.grid.GridCoverage;
055import org.opengis.coverage.grid.GridEnvelope;
056import org.opengis.referencing.FactoryException;
057import org.opengis.referencing.crs.CoordinateReferenceSystem;
058
059import ptolemy.kernel.util.IllegalActionException;
060
061/** A collection of utilities for rasters.
062 * 
063 *  @author Daniel Crawl
064 *  @version $Id: RasterUtilities.java 34712 2018-10-16 03:23:34Z crawl $
065 *  
066 */
067public class RasterUtilities {
068
069    /** This class cannot be instantiated. */
070    private RasterUtilities() { }
071          
072    /** Get a coverage for a raster file.
073     *  @return Returns null if the file is an unknown format.
074     */
075    public static GridCoverage getCoverage(File file, CoordinateReferenceSystem defaultCRS) throws IllegalArgumentException, IOException {
076        AbstractGridCoverage2DReader reader = getReader(file, defaultCRS);
077        if(reader != null) {
078            return reader.read(null);
079        }
080        return null;
081    }
082    
083    /** Get the grid size of a raster.
084     * @return The width and height.
085     */
086    public static Integer[] getGridSize(RasterToken raster) throws IllegalActionException {
087        
088        Integer[] sizes = new Integer[2];
089        
090        File rasterFile = raster.rasterFile();
091        if(rasterFile != null) {
092            throw new IllegalActionException("Getting cell size from file not supported. Use RasterReader to first read the raster.");
093        }
094        
095        GridEnvelope grid = raster.reader().getOriginalGridRange();
096        sizes[0] = grid.getSpan(0);
097        sizes[1] = grid.getSpan(1);
098        
099        return sizes;
100    }
101
102    /** Get a reader for a raster file.
103     *  @return Returns null if file is unknown format.
104     */
105    public static AbstractGridCoverage2DReader getReader(File file, CoordinateReferenceSystem defaultCRS) {
106        AbstractGridFormat format = GridFormatFinder.findFormat(file);
107        if(!(format instanceof UnknownFormat)) {
108
109            AbstractGridCoverage2DReader reader = format.getReader(file);
110            if(defaultCRS == null) {            
111                reader = format.getReader(file);
112            } else {
113                Hints hints = new Hints();
114                hints.put(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, defaultCRS);
115                reader = format.getReader(file, hints);
116            }
117            
118            if(reader != null) {
119                _readers.add(reader);
120            } else {
121                System.err.println("ERROR: null reader for " + file);
122            }
123            return reader;
124        }
125        return null;
126    }
127    
128    /** Close all readers opened by getReader() and getCoverage(). */
129    public static void closeReaders() {
130        
131        synchronized(_readers) {
132            for(AbstractGridCoverage2DReader reader : _readers) {
133                reader.dispose();
134            }   
135            _readers.clear();
136        }        
137    }
138    
139    /** Convert raster file formats using gdal_translate. */
140    public static boolean convertFileFormat(File inputFile, File outputFile,
141        Integer noData) {
142        
143        List<String> args = new ArrayList<String>();
144        args.add("gdal_translate");
145
146        
147        if(noData != null) {
148            args.add("-a_nodata");
149            args.add(noData.toString());
150        }   
151        
152        args.add(inputFile.getAbsolutePath());
153        args.add(outputFile.getAbsolutePath());
154
155        ProcessBuilder builder = new ProcessBuilder(args);
156        builder.redirectErrorStream(true);
157        try {
158            Process process = builder.start();
159            /*
160            try(InputStream stdout = process.getInputStream();
161                    BufferedReader reader = new BufferedReader(new InputStreamReader(stdout));) {
162                    String line;
163                    while((line = reader.readLine()) != null) {
164                        System.out.println(line);
165                    }
166            } finally {
167            */
168            process.waitFor();
169            //}
170        } catch (InterruptedException | IOException e) {
171            System.err.println("WARNING: gdal_translate " +
172                inputFile.getAbsolutePath() + " " + outputFile.getAbsolutePath() +
173                ": " + e.getMessage());
174            return false;
175        }        
176        return true;
177    }
178    
179    /** Get a RasterInfo object for a raster file. This method calls "gdalinfo"
180     *  to find the CRS and bounding box.
181     */
182    public static RasterInfo getInfo(File file) throws IllegalActionException {
183        
184        if(!file.exists()) {
185            throw new IllegalActionException("Raster file does not exist: " + file);
186        }
187        
188        try {
189            ProcessBuilder builder = new ProcessBuilder("gdalinfo", file.getAbsolutePath());
190            builder.redirectErrorStream(true);
191            Process process = builder.start();
192            try(InputStream stdout = process.getInputStream();
193                BufferedReader reader = new BufferedReader(new InputStreamReader(stdout));) {
194                String line;
195                boolean crs = false;
196                boolean upperLeft = false;
197                boolean lowerRight = false;
198                StringBuilder buf = null;
199                double[] minCoords = new double[2];
200                double[] maxCoords = new double[2];
201                RasterInfo info = new RasterInfo();
202                while((line = reader.readLine()) != null) {
203                    
204                    if(line.startsWith("Coordinate System is:")) {
205                        crs = true;
206                        buf = new StringBuilder();
207                    } else if(line.startsWith("Origin = ")) {
208                        crs = false;
209                        if(buf == null) {
210                            System.err.println("WARNING: missing WKT for " + file.getAbsolutePath() +
211                                ". Does .prj file exist?");
212                        } else {
213                            try {
214                                // FIXME buf can be null?
215                                info._crs = CRS.parseWKT(buf.toString());
216                            } catch (FactoryException e) {
217                                System.err.println("WARNING: unable to parse WKT: " + e.getMessage());
218                            }
219                        }
220                    } else if(line.startsWith("GCP[")) {
221                        info._hasGCP = true;
222                    } else if(crs) {
223                        buf.append(line);
224                    } else {
225                        Matcher matcher = UPPER_LEFT_PATTERN.matcher(line);
226                        if(matcher.find()) {
227                            upperLeft = true;
228                            minCoords[0] = Double.parseDouble(matcher.group(1));
229                            maxCoords[1] = Double.parseDouble(matcher.group(2));
230                            continue;
231                        }
232                        matcher = LOWER_RIGHT_PATTERN.matcher(line);
233                        if(matcher.find()) {
234                            lowerRight = true;
235                            maxCoords[0] = Double.parseDouble(matcher.group(1));
236                            minCoords[1] = Double.parseDouble(matcher.group(2));
237                            continue;
238                        }
239                    }
240                }
241                
242                if(upperLeft && lowerRight) {
243                    info._envelope = ReferencedEnvelope.create(
244                            new GeneralEnvelope(minCoords, maxCoords),
245                            info._crs);
246                    //System.out.println("envelope is = " + info._envelope);
247                }
248
249                return info;
250
251            } finally {
252                process.waitFor();
253            }
254            
255        } catch (IOException e) {
256            throw new IllegalActionException("Error executing gdalinfo: " + e.getMessage());
257        } catch (InterruptedException e) {
258            throw new IllegalActionException("Error waiting for gdalinfo to finish: " + e.getMessage());
259        }
260    }
261
262    /** Reproject a raster file.
263     *  @param rasterFile the raster file.
264     *  @param destCRS the CRS to reproject to.
265     */
266    public static File reproject(File rasterFile, CoordinateReferenceSystem destCRS) throws IllegalActionException {
267        
268        // make sure projection is different
269        RasterInfo info = getInfo(rasterFile);
270        if(info.getCRS() == null) {
271            throw new IllegalActionException("Unknown projection for source raster file.");
272        } else if(info.getCRS().equals(destCRS)) {
273            return rasterFile;
274        } else {
275
276            int projectionNum;
277            try {
278                projectionNum = CRS.lookupEpsgCode(destCRS, false);
279            } catch (FactoryException e) {
280                throw new IllegalActionException("Error looking up EPSG code for destination projectin.");
281            }
282            
283            String srcName = FilenameUtils.getBaseName(rasterFile.getAbsolutePath());
284            
285            // FIXME this assumes input raster is arc ascii
286            
287            String destFileNameVRT = rasterFile.getParent() +
288                File.separator +
289                srcName +
290                "-" +
291                projectionNum +
292                ".vrt";
293            
294            ProcessBuilder builder = new ProcessBuilder("gdalwarp",
295                "-t_srs",
296                "EPSG:" + projectionNum,
297                "-of",
298                "VRT",
299                rasterFile.getAbsolutePath(),
300                destFileNameVRT);
301            
302            builder.redirectErrorStream(true);
303            try {
304                Process process = builder.start();
305                /*
306                try {
307                    try(InputStream stdout = process.getInputStream();
308                        BufferedReader reader = new BufferedReader(new InputStreamReader(stdout));) {
309                        String line;
310                        while((line = reader.readLine()) != null) {
311                            System.out.println(line);
312                        }
313                    }
314                } finally {
315                */
316                    if(process.waitFor() != 0) {
317                        throw new IllegalActionException("gdalwarp did not return 0.");
318                    }
319                //}
320            } catch(IOException | InterruptedException e) {
321                throw new IllegalActionException("Error running gdalwarp: " +
322                        e.getMessage());
323            }
324                            
325            String destFileName = rasterFile.getParent() +
326                File.separator +
327                srcName +
328                "-" +
329                projectionNum +
330                "." +
331                FilenameUtils.getExtension(rasterFile.getAbsolutePath());
332            
333            builder = new ProcessBuilder("gdal_translate",
334                "-of",
335                "AAIGrid",
336                destFileNameVRT,
337                destFileName);
338            
339            builder.redirectErrorStream(true);
340            try {
341                Process process = builder.start();
342                /*
343                try {
344                    try(InputStream stdout = process.getInputStream();
345                        BufferedReader reader = new BufferedReader(new InputStreamReader(stdout));) {
346                        String line;
347                        while((line = reader.readLine()) != null) {
348                            System.out.println(line);
349                        }
350                    }
351                } finally {
352                */
353                    if(process.waitFor() != 0) {
354                        throw new IllegalActionException("gdalwarp did not return 0.");
355                    }
356                //}
357            } catch(IOException | InterruptedException e) {
358                throw new IllegalActionException("Error running gdalwarp: " +
359                        e.getMessage());
360            }
361            return new File(destFileName);
362        }              
363    }
364    
365    /** Reproject a raster using gdalwarp to resolve ground control points. */
366    public static void reprojectToResolveGCP(File inputFile, File outputFile) throws IllegalActionException {
367
368        
369        ProcessBuilder builder = new ProcessBuilder("gdalwarp",
370            inputFile.getAbsolutePath(), outputFile.getAbsolutePath());
371        builder.redirectErrorStream(true);
372        try {
373            Process process = builder.start();
374            /*
375            try(InputStream stdout = process.getInputStream();
376                    BufferedReader reader = new BufferedReader(new InputStreamReader(stdout));) {
377                    String line;
378                    while((line = reader.readLine()) != null) {
379                        System.out.println(line);
380                    }
381            } finally {
382            */
383                process.waitFor();
384            //}
385        } catch (IOException | InterruptedException e) {
386            throw new IllegalActionException("Error resolving GCP: " + e.getMessage());
387        }
388        
389    }
390
391
392    /** Utility class that contains information about a raster. */
393    public static class RasterInfo {
394        
395        private RasterInfo() {
396        }
397        
398        public CoordinateReferenceSystem getCRS() {
399            return _crs;
400        }
401        
402        public ReferencedEnvelope getEnvelope() {
403            return _envelope;
404        }
405        
406        public boolean hasGCP() {
407            return _hasGCP;
408        }
409     
410        private CoordinateReferenceSystem _crs;
411        private ReferencedEnvelope _envelope;
412        private boolean _hasGCP;
413    }
414    
415    /** Regex for upper left of bounding box printed by gdalinfo. */
416    private final static Pattern UPPER_LEFT_PATTERN =
417        Pattern.compile("Upper Left\\s+\\(\\s*(\\-?\\d+\\.\\d+),\\s*(\\-?\\d+\\.\\d+)");
418
419    /** Regex for lower right of bounding box printed by gdalinfo. */
420    private final static Pattern LOWER_RIGHT_PATTERN =
421        Pattern.compile("Lower Right\\s+\\(\\s*(\\-?\\d+\\.\\d+),\\s*(\\-?\\d+\\.\\d+)");
422    
423    /** Set of opened readers. */
424    private static Set<AbstractGridCoverage2DReader> _readers = Collections.synchronizedSet(new HashSet<AbstractGridCoverage2DReader>());
425
426}