diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java new file mode 100644 index 0000000..67346d2 --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertSrtmTile.java @@ -0,0 +1,302 @@ +package btools.mapcreator; + +import java.io.*; +import java.util.zip.*; + +public class ConvertSrtmTile +{ + public static int NROWS; + public static int NCOLS; + + public static final short NODATA2 = -32767; // bil-formats nodata + public static final short NODATA = Short.MIN_VALUE; + + static short[] imagePixels; + + public static int[] diffs = new int[100]; + + private static void readBilZip( String filename, int rowOffset, int colOffset, boolean halfCols ) throws Exception + { + int fileRows = 3601; + int fileCols = 3601; + ZipInputStream zis = new ZipInputStream( new BufferedInputStream( new FileInputStream( filename ) ) ); + try + { + for ( ;; ) + { + ZipEntry ze = zis.getNextEntry(); + if ( ze.getName().endsWith( ".bil" ) ) + { + readBilFromStream( zis, rowOffset, colOffset, fileRows, fileCols, halfCols ); + return; + } + } + } + finally + { + zis.close(); + } + } + + private static void readBilFromStream( InputStream is, int rowOffset, int colOffset, int fileRows, int fileCols, boolean halfCols ) + throws Exception + { + DataInputStream dis = new DataInputStream( new BufferedInputStream( is ) ); + for ( int ir = 0; ir < fileRows; ir++ ) + { + int row = rowOffset + ir; + + short lastVal = 0; + boolean fillGap = false; + + for ( int ic = 0; ic < fileCols; ic++ ) + { + int col = colOffset + ic; + short val; + if ( ( ic % 2 ) == 1 && halfCols ) + { + fillGap = true; + } + else + { + int i0 = dis.read(); + int i1 = dis.read(); + + if ( i0 == -1 || i1 == -1 ) + throw new RuntimeException( "unexcepted end of file reading bil entry!" ); + + val = (short) ( ( i1 << 8 ) | i0 ); + + if ( val == NODATA2 ) + { + val = NODATA; + } + + if ( fillGap ) + { + setPixel( row, col - 1, val, lastVal ); + fillGap = false; + } + +if ( row == 18010 ) System.out.print( val + " " ); + + setPixel( row, col, val, val ); + lastVal = val; + } + } + } + } + + private static void setPixel( int row, int col, short val1, short val2 ) + { + if ( row >= 0 && row < NROWS && col >= 0 && col < NCOLS ) + { + if ( val1 != NODATA && val2 != NODATA ) + { + int val = val1 + val2; + if ( val < -32768 || val > 32767 ) + throw new IllegalArgumentException( "val1=" + val1 + " val2=" + val2 ); + + imagePixels[row * NCOLS + col] = (short) ( val ); + } + } + } + + public static void main( String[] args ) throws Exception + { + doConvert( args[0], Integer.parseInt( args[1] ), Integer.parseInt( args[2] ), args[3], null ); + } + + public static void doConvert( String inputDir, int lonDegreeStart, int latDegreeStart, String outputFile, SrtmRaster raster90 ) throws Exception + { + int extraBorder = 10; + int datacells = 0; + int matchingdatacells = 0; + + NROWS = 5 * 3600 + 1 + 2 * extraBorder; + NCOLS = 5 * 3600 + 1 + 2 * extraBorder; + + imagePixels = new short[NROWS * NCOLS]; // 650 MB ! + + // prefill as NODATA + for ( int row = 0; row < NROWS; row++ ) + { + for ( int col = 0; col < NCOLS; col++ ) + { + imagePixels[row * NCOLS + col] = NODATA; + } + } + + for ( int latIdx = -1; latIdx <= 5; latIdx++ ) + { + int latDegree = latDegreeStart + latIdx; + int rowOffset = extraBorder + ( 4 - latIdx ) * 3600; + + for ( int lonIdx = -1; lonIdx <= 5; lonIdx++ ) + { + int lonDegree = lonDegreeStart + lonIdx; + int colOffset = extraBorder + lonIdx * 3600; + + String filename = inputDir + "/" + formatLat( latDegree ) + "_" + formatLon( lonDegree ) + "_1arc_v3_bil.zip"; + File f = new File( filename ); + if ( f.exists() && f.length() > 0 ) + { + System.out.println( "exist: " + filename ); + boolean halfCol = latDegree >= 50 || latDegree < -50; + readBilZip( filename, rowOffset, colOffset, halfCol ); + } + else + { + System.out.println( "none : " + filename ); + } + } + } + + if ( raster90 != null ) + { + for ( int row90 = 0; row90 < 6001; row90++ ) + { + int crow = 3 * row90 + extraBorder; // center row of 3x3 + for ( int col90 = 0; col90 < 6001; col90++ ) + { + int ccol = 3 * col90 + extraBorder; // center col of 3x3 + + short v90 = raster90.eval_array[row90 * 6001 + col90]; + + // evaluate 3x3 area + int sum = 0; + int nodatas = 0; + for ( int row = crow - 1; row <= crow + 1; row++ ) + { + for ( int col = ccol - 1; col <= ccol + 1; col++ ) + { + short v30 = imagePixels[row * NCOLS + col]; + sum += v30; + if ( v30 == NODATA ) + { + nodatas++; + } + } + } + boolean doReplace = nodatas > 0 || v90 == NODATA; + short replaceValue = NODATA; + if ( !doReplace ) + { + datacells++; + int diff = sum - 9 * v90; + replaceValue= (short)(diff / 9); + doReplace = true; + if ( diff < -9 || diff > 9 ) + { +// System.out.println( "replacing due to median missmatch: sum=" + sum + " v90=" + v90 ); + doReplace = true; + } + if ( diff > -50 && diff < 50 ) + { + matchingdatacells++; + diffs[diff+50]++; + } + } + if ( doReplace ) + { + for ( int row = crow - 1; row <= crow + 1; row++ ) + { + for ( int col = ccol - 1; col <= ccol + 1; col++ ) + { +// imagePixels[row * NCOLS + col] = v90; + imagePixels[row * NCOLS + col] = replaceValue; + } + } + } + } + } + } + + SrtmRaster raster = new SrtmRaster(); + raster.nrows = NROWS; + raster.ncols = NCOLS; + raster.noDataValue = NODATA; + raster.cellsize = 1 / 3600.; + raster.xllcorner = lonDegreeStart - ( 0.5 + extraBorder ) * raster.cellsize; + raster.yllcorner = latDegreeStart - ( 0.5 + extraBorder ) * raster.cellsize; + raster.eval_array = imagePixels; + + // encode the raster + OutputStream os = new BufferedOutputStream( new FileOutputStream( outputFile ) ); + new RasterCoder().encodeRaster( raster, os ); + os.close(); + + // decode the raster + InputStream is = new BufferedInputStream( new FileInputStream( outputFile ) ); + SrtmRaster raster2 = new RasterCoder().decodeRaster( is ); + is.close(); + + short[] pix2 = raster2.eval_array; + if ( pix2.length != imagePixels.length ) + throw new RuntimeException( "length mismatch!" ); + + for ( int i = 0; i < pix2.length; i++ ) + { + if ( pix2[i] != imagePixels[i] ) + { + throw new RuntimeException( "content mismatch!" ); + } + } + + for(int i=0; i<100;i++) System.out.println( "diff[" + (i-50) + "] = " + diffs[i] ); + System.out.println( "datacells=" + datacells + " mismatch%=" + 100.*(datacells-matchingdatacells)/datacells ); + // test( raster ); + // raster.calcWeights( 50. ); + // test( raster ); + // 39828330 &lon=3115280&layer=OpenStreetMap + } + + private static void test( SrtmRaster raster ) + { + int lat0 = 39828330; + int lon0 = 3115280; + + for ( int iy = -9; iy <= 9; iy++ ) + { + StringBuilder sb = new StringBuilder(); + for ( int ix = -9; ix <= 9; ix++ ) + { + int lat = lat0 + 90000000 - 100 * iy; + int lon = lon0 + 180000000 + 100 * ix; + int ival = (int) ( raster.getElevation( lon, lat ) / 4. ); + String sval = " " + ival; + sb.append( sval.substring( sval.length() - 4 ) ); + } + System.out.println( sb ); + System.out.println(); + } + } + + private static String formatLon( int lon ) + { + if ( lon >= 180 ) + lon -= 180; // TODO: w180 oder E180 ? + + String s = "e"; + if ( lon < 0 ) + { + lon = -lon; + s = "w"; + } + String n = "000" + lon; + return s + n.substring( n.length() - 3 ); + } + + private static String formatLat( int lat ) + { + String s = "n"; + if ( lat < 0 ) + { + lat = -lat; + s = "s"; + } + String n = "00" + lat; + return s + n.substring( n.length() - 2 ); + } + +} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java new file mode 100644 index 0000000..c2ba42b --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/ConvertUrlList.java @@ -0,0 +1,65 @@ +package btools.mapcreator; + +import java.io.BufferedReader; +import java.io.File; +import java.io.FileReader; + +public class ConvertUrlList +{ + public static final short NODATA = -32767; + + public static void main( String[] args ) throws Exception + { + BufferedReader br = new BufferedReader( new FileReader( args[0] ) ); + + for ( ;; ) + { + String line = br.readLine(); + if ( line == null ) + { + break; + } + int idx1 = line.indexOf( "srtm_" ); + if ( idx1 < 0 ) + { + continue; + } + + String filename90 = line.substring( idx1 ); + String filename30 = filename90.substring( 0, filename90.length() - 3 ) + "bef"; + + if ( new File( filename30 ).exists() ) + { + continue; + } + + // int srtmLonIdx = (ilon+5000000)/5000000; -> ilon = (srtmLonIdx-1)*5 + // int srtmLatIdx = (154999999-ilat)/5000000; -> ilat = 155 - srtmLatIdx*5 + + int srtmLonIdx = Integer.parseInt( filename90.substring( 5, 7 ).toLowerCase() ); + int srtmLatIdx = Integer.parseInt( filename90.substring( 8, 10 ).toLowerCase() ); + + int ilon_base = ( srtmLonIdx - 1 ) * 5 - 180; + int ilat_base = 150 - srtmLatIdx * 5 - 90; + + if ( ilon_base < 5 || ilon_base > 10 || ilat_base < 45 || ilat_base > 50 ) + { + continue; + } + + SrtmRaster raster90 = null; + + File file90 = new File( new File( args[1] ), filename90 ); + if ( file90.exists() ) + { + System.out.println( "reading " + file90 ); + raster90 = new SrtmData( file90 ).getRaster(); + } + + ConvertSrtmTile.doConvert( args[2], ilon_base, ilat_base, filename30, raster90 ); + } + br.close(); + } + + +} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java b/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java index 205a5b9..b45aec1 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/PosUnifier.java @@ -1,9 +1,11 @@ package btools.mapcreator; +import java.io.BufferedInputStream; import java.io.DataInputStream; -import java.io.DataOutputStream; import java.io.EOFException; import java.io.File; +import java.io.FileInputStream; +import java.io.InputStream; import java.util.HashMap; import btools.util.CompactLongSet; @@ -13,10 +15,8 @@ import btools.util.FrozenLongSet; /** * PosUnifier does 3 steps in map-processing: * - * - unify positions - * - add srtm elevation data - * - make a bordernodes file containing net data - * from the bordernids-file just containing ids + * - unify positions - add srtm elevation data - make a bordernodes file + * containing net data from the bordernids-file just containing ids * * @author ab */ @@ -27,21 +27,20 @@ public class PosUnifier extends MapCreatorBase private File nodeTilesOut; private CompactLongSet positionSet; - private HashMap srtmmap ; + private HashMap srtmmap; private int lastStrmLonIdx; private int lastStrmLatIdx; - private SrtmData lastSrtmData; + private SrtmRaster lastSrtmRaster; private String srtmdir; private CompactLongSet borderNids; - - public static void main(String[] args) throws Exception + public static void main( String[] args ) throws Exception { - System.out.println("*** PosUnifier: Unify position values and enhance elevation"); - if (args.length != 5) + System.out.println( "*** PosUnifier: Unify position values and enhance elevation" ); + if ( args.length != 5 ) { - System.out.println("usage: java PosUnifier " ); + System.out.println( "usage: java PosUnifier " ); return; } new PosUnifier().process( new File( args[0] ), new File( args[1] ), new File( args[2] ), new File( args[3] ), args[4] ); @@ -57,13 +56,14 @@ public class PosUnifier extends MapCreatorBase borderNids = new CompactLongSet(); try { - for(;;) + for ( ;; ) { long nid = readId( dis ); - if ( !borderNids.contains( nid ) ) borderNids.fastAdd( nid ); + if ( !borderNids.contains( nid ) ) + borderNids.fastAdd( nid ); } } - catch( EOFException eof ) + catch (EOFException eof) { dis.close(); } @@ -88,8 +88,8 @@ public class PosUnifier extends MapCreatorBase @Override public void nextNode( NodeData n ) throws Exception { - SrtmData srtm = srtmForNode( n.ilon, n.ilat ); - n.selev = srtm == null ? Short.MIN_VALUE : srtm.getElevation( n.ilon, n.ilat); + SrtmRaster srtm = srtmForNode( n.ilon, n.ilat ); + n.selev = srtm == null ? Short.MIN_VALUE : srtm.getElevation( n.ilon, n.ilat ); findUniquePos( n ); @@ -113,13 +113,13 @@ public class PosUnifier extends MapCreatorBase int londelta = lonmod < 500000 ? 1 : -1; int latmod = n.ilat % 1000000; int latdelta = latmod < 500000 ? 1 : -1; - for(int latsteps = 0; latsteps < 100; latsteps++) + for ( int latsteps = 0; latsteps < 100; latsteps++ ) { - for(int lonsteps = 0; lonsteps <= latsteps; lonsteps++) + for ( int lonsteps = 0; lonsteps <= latsteps; lonsteps++ ) { - int lon = n.ilon + lonsteps*londelta; - int lat = n.ilat + latsteps*latdelta; - long pid = ((long)lon)<<32 | lat; // id from position + int lon = n.ilon + lonsteps * londelta; + int lat = n.ilat + latsteps * latdelta; + long pid = ( (long) lon ) << 32 | lat; // id from position if ( !positionSet.contains( pid ) ) { positionSet.fastAdd( pid ); @@ -132,16 +132,14 @@ public class PosUnifier extends MapCreatorBase System.out.println( "*** WARNING: cannot unify position for: " + n.ilon + " " + n.ilat ); } - /** - * get the srtm data set for a position - * srtm coords are srtm__ - * where srtmLon = 180 + lon, srtmLat = 60 - lat + * get the srtm data set for a position srtm coords are + * srtm__ where srtmLon = 180 + lon, srtmLat = 60 - lat */ - private SrtmData srtmForNode( int ilon, int ilat ) throws Exception + private SrtmRaster srtmForNode( int ilon, int ilat ) throws Exception { - int srtmLonIdx = (ilon+5000000)/5000000; - int srtmLatIdx = (154999999-ilat)/5000000; + int srtmLonIdx = ( ilon + 5000000 ) / 5000000; + int srtmLatIdx = ( 154999999 - ilat ) / 5000000; if ( srtmLatIdx < 1 || srtmLatIdx > 24 || srtmLonIdx < 1 || srtmLonIdx > 72 ) { @@ -149,45 +147,63 @@ public class PosUnifier extends MapCreatorBase } if ( srtmLonIdx == lastStrmLonIdx && srtmLatIdx == lastStrmLatIdx ) { - return lastSrtmData; + return lastSrtmRaster; } lastStrmLonIdx = srtmLonIdx; lastStrmLatIdx = srtmLatIdx; StringBuilder sb = new StringBuilder( 16 ); sb.append( "srtm_" ); - sb.append( (char)('0' + srtmLonIdx/10 ) ).append( (char)('0' + srtmLonIdx%10 ) ).append( '_' ); - sb.append( (char)('0' + srtmLatIdx/10 ) ).append( (char)('0' + srtmLatIdx%10 ) ).append( ".zip" ); + sb.append( (char) ( '0' + srtmLonIdx / 10 ) ).append( (char) ( '0' + srtmLonIdx % 10 ) ).append( '_' ); + sb.append( (char) ( '0' + srtmLatIdx / 10 ) ).append( (char) ( '0' + srtmLatIdx % 10 ) ); String filename = sb.toString(); - - lastSrtmData = srtmmap.get( filename ); - if ( lastSrtmData == null && !srtmmap.containsKey( filename ) ) + lastSrtmRaster = srtmmap.get( filename ); + if ( lastSrtmRaster == null && !srtmmap.containsKey( filename ) ) { - File f = new File( new File( srtmdir ), filename ); + File f = new File( new File( srtmdir ), filename + ".bef" ); + System.out.println( "checking: " + f + " ilon=" + ilon + " ilat=" + ilat ); + if ( f.exists() ) + { + System.out.println( "*** reading: " + f ); + try + { + InputStream isc = new BufferedInputStream( new FileInputStream( f ) ); + lastSrtmRaster = new RasterCoder().decodeRaster( isc ); + isc.close(); + } + catch (Exception e) + { + System.out.println( "**** ERROR reading " + f + " ****" ); + } + srtmmap.put( filename, lastSrtmRaster ); + return lastSrtmRaster; + } + + f = new File( new File( srtmdir ), filename + ".zip" ); System.out.println( "reading: " + f + " ilon=" + ilon + " ilat=" + ilat ); if ( f.exists() ) { - try - { - lastSrtmData = new SrtmData( f ); - } - catch( Exception e ) - { - System.out.println( "**** ERROR reading " + f + " ****" ); - } + try + { + lastSrtmRaster = new SrtmData( f ).getRaster(); + } + catch (Exception e) + { + System.out.println( "**** ERROR reading " + f + " ****" ); + } } - srtmmap.put( filename, lastSrtmData ); + srtmmap.put( filename, lastSrtmRaster ); } - return lastSrtmData; + return lastSrtmRaster; } private void resetSrtm() { - srtmmap = new HashMap(); + srtmmap = new HashMap(); lastStrmLonIdx = -1; lastStrmLatIdx = -1; - lastSrtmData = null; + lastSrtmRaster = null; } } diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java new file mode 100644 index 0000000..4471612 --- /dev/null +++ b/brouter-map-creator/src/main/java/btools/mapcreator/RasterCoder.java @@ -0,0 +1,88 @@ +package btools.mapcreator; + +import java.io.*; +import btools.util.*; + +// +// Encode/decode a raster +// + +public class RasterCoder +{ + public void encodeRaster(SrtmRaster raster, OutputStream os) throws IOException + { + DataOutputStream dos = new DataOutputStream(os); + + long t0 = System.currentTimeMillis(); + + dos.writeInt(raster.ncols); + dos.writeInt(raster.nrows); + dos.writeDouble(raster.xllcorner); + dos.writeDouble(raster.yllcorner); + dos.writeDouble(raster.cellsize); + dos.writeShort(raster.noDataValue); + + _encodeRaster(raster, os); + long t1 = System.currentTimeMillis(); + + System.out.println("finished encoding in " + (t1 - t0) + " ms"); + } + + public SrtmRaster decodeRaster(InputStream is) throws IOException + { + DataInputStream dis = new DataInputStream(is); + + long t0 = System.currentTimeMillis(); + + SrtmRaster raster = new SrtmRaster(); + raster.ncols = dis.readInt(); + raster.nrows = dis.readInt(); + raster.xllcorner = dis.readDouble(); + raster.yllcorner = dis.readDouble(); + raster.cellsize = dis.readDouble(); + raster.noDataValue = dis.readShort(); + raster.eval_array = new short[raster.ncols * raster.nrows]; + + _decodeRaster(raster, is); + + raster.usingWeights = true; + + long t1 = System.currentTimeMillis(); + System.out.println("finished decoding in " + (t1 - t0) + " ms"); + return raster; + } + + + private void _encodeRaster(SrtmRaster raster, OutputStream os) throws IOException + { + MixCoderDataOutputStream mco = new MixCoderDataOutputStream(os); + int nrows = raster.nrows; + int ncols = raster.ncols; + short[] pixels = raster.eval_array; + + for (int row = 0; row < nrows; row++) + { + for (int col = 0; col < ncols; col++) + { + mco.writeMixed(pixels[row * ncols + col]); + } + } + mco.flush(); + } + + private void _decodeRaster(SrtmRaster raster, InputStream is) throws IOException + { + MixCoderDataInputStream mci = new MixCoderDataInputStream(is); + int nrows = raster.nrows; + int ncols = raster.ncols; + short[] pixels = raster.eval_array; + + for (int row = 0; row < nrows; row++) + { + for (int col = 0; col < ncols; col++) + { + pixels[row * ncols + col] = (short) mci.readMixed(); + } + } + } +} diff --git a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java index abcc6fc..ea42f96 100644 --- a/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java +++ b/brouter-map-creator/src/main/java/btools/mapcreator/SrtmData.java @@ -1,3 +1,5 @@ +package btools.mapcreator; + /** * This is a wrapper for a 5*5 degree srtm file in ascii/zip-format * @@ -7,8 +9,8 @@ * * @author ab */ -package btools.mapcreator; +import java.io.BufferedInputStream; import java.io.BufferedReader; import java.io.File; import java.io.FileInputStream; @@ -18,63 +20,18 @@ import java.util.StringTokenizer; import java.util.zip.ZipEntry; import java.util.zip.ZipInputStream; - public class SrtmData { - public int ncols; - public int nrows; - public double xllcorner; - public double yllcorner; - public double cellsize; - public short[] eval_array; - - private double minlon; - private double minlat; - - public void init() - { - minlon = xllcorner; - minlat = yllcorner; - } - - private boolean missingData = false; - - public short getElevation( int ilon, int ilat ) - { - double lon = ilon / 1000000. - 180.; - double lat = ilat / 1000000. - 90.; - - double dcol = (lon - minlon)/cellsize -0.5; - double drow = (lat - minlat)/cellsize -0.5; - int row = (int)drow; - int col = (int)dcol; - if ( col < 0 ) col = 0; - if ( col >= ncols-1 ) col = ncols - 2; - if ( row < 0 ) row = 0; - if ( row >= nrows-1 ) row = nrows - 2; - double wrow = drow-row; - double wcol = dcol-col; - missingData = false; - double eval = (1.-wrow)*(1.-wcol)*get(row ,col ) - + ( wrow)*(1.-wcol)*get(row+1,col ) - + (1.-wrow)*( wcol)*get(row ,col+1) - + ( wrow)*( wcol)*get(row+1,col+1); - return missingData ? Short.MIN_VALUE : (short)(eval); - } - - private short get( int r, int c ) - { - short e = eval_array[r*ncols + c ]; - if ( e == Short.MIN_VALUE ) missingData = true; - return e; - } + private SrtmRaster raster; public SrtmData( File file ) throws Exception { - ZipInputStream zis = new ZipInputStream( new FileInputStream( file ) ); + raster = new SrtmRaster(); + + ZipInputStream zis = new ZipInputStream( new BufferedInputStream( new FileInputStream( file ) ) ); try { - for(;;) + for ( ;; ) { ZipEntry ze = zis.getNextEntry(); if ( ze.getName().endsWith( ".asc" ) ) @@ -90,6 +47,11 @@ public class SrtmData } } + public SrtmRaster getRaster() + { + return raster; + } + private String secondToken( String s ) { StringTokenizer tk = new StringTokenizer( s, " " ); @@ -99,23 +61,29 @@ public class SrtmData public void readFromStream( InputStream is ) throws Exception { - BufferedReader br = new BufferedReader(new InputStreamReader(is)); + BufferedReader br = new BufferedReader( new InputStreamReader( is ) ); int linenr = 0; - for(;;) + for ( ;; ) { linenr++; if ( linenr <= 6 ) { String line = br.readLine(); - if ( linenr == 1 ) ncols = Integer.parseInt( secondToken( line ) ); - else if ( linenr == 2 ) nrows = Integer.parseInt( secondToken( line ) ); - else if ( linenr == 3 ) xllcorner = Double.parseDouble( secondToken( line ) ); - else if ( linenr == 4 ) yllcorner = Double.parseDouble( secondToken( line ) ); - else if ( linenr == 5 ) cellsize = Double.parseDouble( secondToken( line ) ); + if ( linenr == 1 ) + raster.ncols = Integer.parseInt( secondToken( line ) ); + else if ( linenr == 2 ) + raster.nrows = Integer.parseInt( secondToken( line ) ); + else if ( linenr == 3 ) + raster.xllcorner = Double.parseDouble( secondToken( line ) ); + else if ( linenr == 4 ) + raster.yllcorner = Double.parseDouble( secondToken( line ) ); + else if ( linenr == 5 ) + raster.cellsize = Double.parseDouble( secondToken( line ) ); else if ( linenr == 6 ) { - // nodata_value ignored, assumed something << 0 - eval_array = new short[ncols * nrows]; + // nodata ignored here ( < -250 assumed nodata... ) + // raster.noDataValue = Short.parseShort( secondToken( line ) ); + raster.eval_array = new short[raster.ncols * raster.nrows]; } } else @@ -124,17 +92,19 @@ public class SrtmData int col = 0; int n = 0; boolean negative = false; - for(;;) + for ( ;; ) { int c = br.read(); - if ( c < 0 ) break; + if ( c < 0 ) + break; if ( c == ' ' ) { - if ( negative ) n = -n; - short val = n < -250 ? Short.MIN_VALUE : (short)(n*4); + if ( negative ) + n = -n; + short val = n < -250 ? Short.MIN_VALUE : (short) (n*2); - eval_array[ (nrows-1-row)*ncols + col ] = val; - if (++col == ncols ) + raster.eval_array[row * raster.ncols + col] = val; + if ( ++col == raster.ncols ) { col = 0; ++row; @@ -144,7 +114,7 @@ public class SrtmData } else if ( c >= '0' && c <= '9' ) { - n = 10*n + (c-'0'); + n = 10 * n + ( c - '0' ); } else if ( c == '-' ) { @@ -154,27 +124,6 @@ public class SrtmData break; } } - init(); br.close(); } - - private void test() - { - int[] ca = new int[]{ 50477121, 8051915, // 181 - 50477742, 8047408, // 154 - 50477189, 8047308, // 159 - }; - for( int i=0; i 0 ) + { + filterDiscRadius = Integer.parseInt( sRadius ); + System.out.println( "using filterDiscRadius = " + filterDiscRadius ); + } + String sFraction = System.getProperty( "filterCenterFraction" ); + if ( sFraction != null && sFraction.length() > 0 ) + { + filterCenterFraction = Integer.parseInt( sFraction ) / 100.; + System.out.println( "using filterCenterFraction = " + filterCenterFraction ); + } + } + + + // calculate interpolation weights from the overlap of a probe disc of given radius at given latitude + // ( latIndex = 0 -> 0 deg, latIndex = 16 -> 80 degree) + + private static Weights[][] getWeights( int latIndex ) + { + int idx = latIndex < 16 ? latIndex : 16; + + Weights[][] res = allShiftWeights[idx]; + if ( res == null ) + { + res = calcWeights( idx ); + allShiftWeights[idx] = res; + } + return res; + } + + private static Weights[][] calcWeights( int latIndex ) + { + double coslat = Math.cos( latIndex * 5. / 57.3 ); + + // radius in pixel units + double ry = filterDiscRadius; + double rx = ry / coslat; + + // gridsize is 2*radius + 1 cell + int nx = ((int)rx) *2 + 3; + int ny = ((int)ry) *2 + 3; + + System.out.println( "nx="+ nx + " ny=" + ny ); + + int mx = nx / 2; // mean pixels + int my = ny / 2; + + // create a matrix for the relative intergrid-position + + Weights[][] shiftWeights = new Weights[gridSteps+1][]; + + // loop the intergrid-position + for( int gx=0; gx<=gridSteps; gx++ ) + { + shiftWeights[gx] = new Weights[gridSteps+1]; + double x0 = mx + ( (double)gx ) / gridSteps; + + for( int gy=0; gy<=gridSteps; gy++ ) + { + double y0 = my + ( (double)gy ) / gridSteps; + + // create the weight-matrix + Weights weights = new Weights( nx, ny ); + shiftWeights[gx][gy] = weights; + + double sampleStep = 0.001; + + for( double x = -1. + sampleStep/2.; x < 1.; x += sampleStep ) + { + double mx2 = 1. - x*x; + + int x_idx = (int)(x0 + x*rx); + + for( double y = -1. + sampleStep/2.; y < 1.; y += sampleStep ) + { + if ( y*y > mx2 ) + { + continue; + } + // we are in the ellipse, see what pixel we are on + int y_idx = (int)(y0 + y*ry); + weights.inc( x_idx, y_idx ); + } + } + weights.normalize( true ); + } + } + return shiftWeights; + } + +}