srtm-update
This commit is contained in:
Arndt
parent
d277db0bd1
commit
19fd3112bf
6 changed files with 845 additions and 136 deletions
|
|
@ -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 );
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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();
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
|
@ -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<String,SrtmData> srtmmap ;
|
||||
private HashMap<String, SrtmRaster> 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 <node-tiles-in> <node-tiles-out> <bordernids-in> <bordernodes-out> <strm-data-dir>" );
|
||||
System.out.println( "usage: java PosUnifier <node-tiles-in> <node-tiles-out> <bordernids-in> <bordernodes-out> <strm-data-dir>" );
|
||||
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_<srtmLon>_<srtmLat>
|
||||
* where srtmLon = 180 + lon, srtmLat = 60 - lat
|
||||
* get the srtm data set for a position srtm coords are
|
||||
* srtm_<srtmLon>_<srtmLat> 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<String,SrtmData>();
|
||||
srtmmap = new HashMap<String, SrtmRaster>();
|
||||
lastStrmLonIdx = -1;
|
||||
lastStrmLatIdx = -1;
|
||||
lastSrtmData = null;
|
||||
lastSrtmRaster = null;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -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<ca.length; i+=2 )
|
||||
{
|
||||
int lat=ca[i] + 90000000;
|
||||
int lon=ca[i+1] + 180000000;
|
||||
System.err.println( "lat=" + lat + " lon=" + lon + " elev=" + getElevation( lon, lat )/4. );
|
||||
}
|
||||
}
|
||||
|
||||
public static void main( String[] args ) throws Exception
|
||||
{
|
||||
SrtmData data = new SrtmData( new File( args[0] ) );
|
||||
data.test();
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,289 @@
|
|||
package btools.mapcreator;
|
||||
|
||||
import btools.util.ReducedMedianFilter;
|
||||
|
||||
/**
|
||||
* Container for a srtm-raster + it's meta-data
|
||||
*
|
||||
* @author ab
|
||||
*/
|
||||
public class SrtmRaster
|
||||
{
|
||||
public int ncols;
|
||||
public int nrows;
|
||||
public double xllcorner;
|
||||
public double yllcorner;
|
||||
public double cellsize;
|
||||
public short[] eval_array;
|
||||
public short noDataValue;
|
||||
|
||||
public boolean usingWeights = false;
|
||||
|
||||
private boolean missingData = false;
|
||||
|
||||
public short getElevation( int ilon, int ilat )
|
||||
{
|
||||
double lon = ilon / 1000000. - 180.;
|
||||
double lat = ilat / 1000000. - 90.;
|
||||
|
||||
if ( usingWeights )
|
||||
{
|
||||
return getElevationFromShiftWeights( lon, lat );
|
||||
}
|
||||
|
||||
// no weights calculated, use 2d linear interpolation
|
||||
double dcol = (lon - xllcorner)/cellsize -0.5;
|
||||
double drow = (lat - yllcorner)/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;
|
||||
|
||||
// System.out.println( "wrow=" + wrow + " wcol=" + wcol + " row=" + row + " col=" + col );
|
||||
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);
|
||||
// System.out.println( "eval=" + eval );
|
||||
return missingData ? Short.MIN_VALUE : (short)(eval*2);
|
||||
}
|
||||
|
||||
private short get( int r, int c )
|
||||
{
|
||||
short e = eval_array[ (nrows-1-r)*ncols + c ];
|
||||
if ( e == Short.MIN_VALUE ) missingData = true;
|
||||
return e;
|
||||
}
|
||||
|
||||
private short getElevationFromShiftWeights( double lon, double lat )
|
||||
{
|
||||
// calc lat-idx and -weight
|
||||
double alat = lat < 0. ? - lat : lat;
|
||||
alat /= 5.;
|
||||
int latIdx = (int)alat;
|
||||
double wlat = alat - latIdx;
|
||||
|
||||
double dcol = (lon - xllcorner)/cellsize;
|
||||
double drow = (lat - yllcorner)/cellsize;
|
||||
int row = (int)drow;
|
||||
int col = (int)dcol;
|
||||
if ( col < 9 || col >= ncols-9 || row < 9 || row >= nrows-9 )
|
||||
{
|
||||
throw new IllegalArgumentException( "out of boundary range: col=" + col + " row=" + row );
|
||||
}
|
||||
double dgx = (dcol-col)*gridSteps;
|
||||
double dgy = (drow-row)*gridSteps;
|
||||
|
||||
// System.out.println( "wrow=" + wrow + " wcol=" + wcol + " row=" + row + " col=" + col );
|
||||
|
||||
int gx = (int)(dgx);
|
||||
int gy = (int)(dgy);
|
||||
|
||||
double wx = dgx-gx;
|
||||
double wy = dgy-gy;
|
||||
|
||||
double w00 = (1.-wx)*(1.-wy);
|
||||
double w01 = (1.-wx)*( wy);
|
||||
double w10 = ( wx)*(1.-wy);
|
||||
double w11 = ( wx)*( wy);
|
||||
|
||||
Weights[][] w0 = getWeights( latIdx );
|
||||
Weights[][] w1 = getWeights( latIdx+1 );
|
||||
|
||||
missingData = false;
|
||||
|
||||
double m0 = w00*getElevation( w0[gx ][gy ], row, col )
|
||||
+ w01*getElevation( w0[gx ][gy+1], row, col )
|
||||
+ w10*getElevation( w0[gx+1][gy ], row, col )
|
||||
+ w11*getElevation( w0[gx+1][gy+1], row, col );
|
||||
double m1 = w00*getElevation( w1[gx ][gy ], row, col )
|
||||
+ w01*getElevation( w1[gx ][gy+1], row, col )
|
||||
+ w10*getElevation( w1[gx+1][gy ], row, col )
|
||||
+ w11*getElevation( w1[gx+1][gy+1], row, col );
|
||||
|
||||
if ( missingData ) return Short.MIN_VALUE;
|
||||
double m = (1.-wlat) * m0 + wlat * m1;
|
||||
return (short)(m * 2);
|
||||
}
|
||||
|
||||
private double getElevation( Weights w, int row, int col )
|
||||
{
|
||||
if ( missingData )
|
||||
{
|
||||
return 0.;
|
||||
}
|
||||
int nx = w.nx;
|
||||
int ny = w.ny;
|
||||
int mx = nx / 2; // mean pixels
|
||||
int my = ny / 2;
|
||||
|
||||
System.out.println( "nx="+ nx + " ny=" + ny );
|
||||
|
||||
ReducedMedianFilter rmf = new ReducedMedianFilter( nx * ny );
|
||||
|
||||
for( int ix = 0; ix < nx; ix ++ )
|
||||
{
|
||||
for( int iy = 0; iy < ny; iy ++ )
|
||||
{
|
||||
short val = get( row + iy - my, col + ix - mx );
|
||||
rmf.addSample( w.getWeight( ix, iy ), val );
|
||||
}
|
||||
}
|
||||
return missingData ? rmf.calcEdgeReducedMedian( filterCenterFraction ) : 0.;
|
||||
}
|
||||
|
||||
|
||||
private static class Weights
|
||||
{
|
||||
int nx;
|
||||
int ny;
|
||||
double[] weights;
|
||||
long total = 0;
|
||||
|
||||
Weights( int nx, int ny )
|
||||
{
|
||||
this.nx = nx;
|
||||
this.ny = ny;
|
||||
weights = new double[nx*ny];
|
||||
}
|
||||
|
||||
void inc( int ix, int iy )
|
||||
{
|
||||
weights[ iy*nx + ix ] += 1.;
|
||||
total++;
|
||||
}
|
||||
|
||||
void normalize( boolean verbose )
|
||||
{
|
||||
for( int iy =0; iy < ny; iy++ )
|
||||
{
|
||||
StringBuilder sb = verbose ? new StringBuilder() : null;
|
||||
for( int ix =0; ix < nx; ix++ )
|
||||
{
|
||||
weights[ iy*nx + ix ] /= total;
|
||||
if ( sb != null )
|
||||
{
|
||||
int iweight = (int)(1000*weights[ iy*nx + ix ] + 0.5 );
|
||||
String sval = " " + iweight;
|
||||
sb.append( sval.substring( sval.length() - 4 ) );
|
||||
}
|
||||
}
|
||||
if ( sb != null )
|
||||
{
|
||||
System.out.println( sb );
|
||||
System.out.println();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double getWeight( int ix, int iy )
|
||||
{
|
||||
return weights[ iy*nx + ix ];
|
||||
}
|
||||
}
|
||||
|
||||
private static int gridSteps = 10;
|
||||
private static Weights[][][] allShiftWeights = new Weights[17][][];
|
||||
|
||||
private static double filterCenterFraction = 0.5;
|
||||
private static double filterDiscRadius = 2.9; // in pixels
|
||||
|
||||
static
|
||||
{
|
||||
String sRadius = System.getProperty( "filterDiscRadius" );
|
||||
if ( sRadius != null && sRadius.length() > 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;
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Reference in a new issue