crotwell / seisFile

A library for reading and writing seismic file formats in java.
GNU Lesser General Public License v3.0
28 stars 20 forks source link

How to fill gaps? #33

Closed andreabono closed 8 months ago

andreabono commented 8 months ago

I wrote this routine to read miniseed data from an Earthworm waveserver into a local miniseed file.

public String  readData(String station , String channel,  String network, 
            String location, Date begin, Date end, String output_path) throws IOException{

        String outputFile = output_path +station+"."+channel+".mseed";
        DataOutputStream out = null;
        try {     
            query = reader.createQuery(network, station, location, channel, begin.toInstant(), end.toInstant());

            List<DataRecord> dati = reader.read(query);
            if (dati!=null){      
                ListIterator drIter = dati.listIterator();

                if (!drIter.hasNext()) {
                    App.logger.info("No data for " + query);
                    return null;
                } else {
                    out = new DataOutputStream(new BufferedOutputStream(new FileOutputStream(outputFile)));

                    DataRecord dr;
                    while (drIter.hasNext()) {
                        dr = (DataRecord)drIter.next();

                        dr.write(out);
                    }
                    dr=null;

                    return outputFile;
                }
            } else 
                return null;

        } catch (Exception ex) {
            ex.printStackTrace();
            Logger.getLogger("org##### ").log(java.util.logging.Level.SEVERE, ex.getMessage());
            return null;
        } finally {
            try {
                out.flush();
                out.close();
            } catch (Exception ex) {}
        }
    } 

But now I have troubles with waves with gaps. How to fill such gaps with zeroes (or better last good sample)?
I guess the code to be corrected is this one:

while (drIter.hasNext()) {
                        dr = (DataRecord)drIter.next();

                        dr.write(out);
                    }

but I cant find any usefoul solution.

andreabono commented 8 months ago

I'm trying to upgrade my code like this:

public String readData(String station, String channel, String network, String location, Date begin, Date end, String output_path) throws IOException {
            String outputFile = output_path + station + "." + channel + ".mseed";
            DataOutputStream out = null;
            try {
                //reader.setVerbose(true);
                query = reader.createQuery(network, station, location, channel, begin.toInstant(), end.toInstant());

                List<DataRecord> dati = reader.read(query);
                if (dati != null) {
                    ListIterator<DataRecord> drIter = dati.listIterator();
                    if (!drIter.hasNext()) {
                        App.logger.info("No data for " + query);
                        return null;
                    } else {

                        out = new DataOutputStream(new BufferedOutputStream(new FileOutputStream(outputFile)));

                        DataRecord prevRecord = null;
                        int packetduration;
                        while (drIter.hasNext()) {
                            DataRecord currentRecord = drIter.next();
                            // Create aDateTimeFormatter 
                            DateTimeFormatter formatter = DateTimeFormatter.ofPattern("yyyy,DDD,HH:mm:ss.SSSS");
                            // Convert the sing into an Instant
                            Instant currentEndTime = LocalDateTime.parse(currentRecord.getEndTime(), formatter).toInstant(ZoneOffset.UTC);

                            if (prevRecord != null) {
                                packetduration = currentRecord.getHeader().getNumSamples() / Math.abs(currentRecord.getHeader().getSampleRateFactor()/currentRecord.getHeader().getSampleRateMultiplier());

                                Instant prevEndTime = LocalDateTime.parse(prevRecord.getEndTime(), formatter).toInstant(ZoneOffset.UTC);

                                if (ChronoUnit.SECONDS.between(currentEndTime, prevEndTime.plusSeconds(packetduration)) < 0) {
                                    long differenza =ChronoUnit.SECONDS.between(currentEndTime, prevEndTime.plusSeconds(packetduration));
                                    // Calcola quanti campioni mancano tra i due record
                                    long missingSamples = Math.abs(differenza) * Math.abs(currentRecord.getHeader().getSampleRateFactor()/currentRecord.getHeader().getSampleRateMultiplier());

                                 //******************************
                                 // WHAT NOW??
                                // How do I add some zeroes to out?
                                //******************************
..........
                                }
                            }

                            currentRecord.write(out);
                            prevRecord = currentRecord;
                        }
                        return outputFile;
                    }
                } else
                    return null;

            } catch (Exception ex) {
                ex.printStackTrace();
                Logger.getLogger("org.ingv.pfx ").log(java.util.logging.Level.SEVERE, ex.getMessage());
                return null;
            } finally {
                try {
                    if (out != null) {
                        out.flush();
                        out.close();
                    }
                } catch (Exception ex) {
                    ex.printStackTrace();
                }
            }
        }
crotwell commented 8 months ago

Gap filling can be very tricky, especially if the gap is not exactly an integral number of samples in length. I guess you could create new DataRecords and insert fake data, but that is generally not considered a good idea because then the reading application can't tell what data is real and what is gap filling. I would recommend that any gap filling operations be done immediately prior to processing in the application instead of in miniseed. Then you are just dealing with arrays of integers in memory, where it is much easier to fix issues like this.

andreabono commented 8 months ago

Ok, I'll to work in the following parts of my client... Thanks!!

andreabono commented 8 months ago

I solved the issue changing my "data extraction" routine that's used while reading data from a file on disk:

 public float[] extract(List<DataRecord> drList) throws UnsupportedCompressionType, CodecException,
            SeedFormatException {
//
        ArrayList<Float> realsamples = new ArrayList<>();

        try {
            Float lgs=Float.valueOf("0");
            DataRecord prevRecord = null;
            int packetduration;
            DateTimeFormatter formatter = DateTimeFormatter.ofPattern("yyyy,DDD,HH:mm:ss.SSSS");
            for (DataRecord currentRecord : drList) {
                Instant currentEndTime = LocalDateTime.parse(currentRecord.getEndTime(), formatter).toInstant(ZoneOffset.UTC);
                if (prevRecord != null) {
                    packetduration = currentRecord.getHeader().getNumSamples() / Math.abs(currentRecord.getHeader().getSampleRateFactor()/currentRecord.getHeader().getSampleRateMultiplier());
                    Instant prevEndTime = LocalDateTime.parse(prevRecord.getEndTime(), formatter).toInstant(ZoneOffset.UTC);
                    if (ChronoUnit.SECONDS.between(currentEndTime, prevEndTime.plusSeconds(packetduration)) < 0) {
                        //-----------------------------------
                        // There's a GAP in the time series
                        //-----------------------------------
                        long differenza =ChronoUnit.SECONDS.between(currentEndTime, prevEndTime.plusSeconds(packetduration));
                        long missingSamples = Math.abs(differenza) * Math.abs(currentRecord.getHeader().getSampleRateFactor()/currentRecord.getHeader().getSampleRateMultiplier());
//                      System.out.println(">>>>>>>>>>>>>>>>>>>> Campioni mancanti " + missingSamples + " su " + 
//                          currentRecord.getHeader().getStationIdentifier() + " " + currentRecord.getHeader().getChannelIdentifier());
                        for (int i =0; i< missingSamples; i++) realsamples.add(Float.valueOf(lgs));
                    }
                }

                //--------------------------------------------
                // Reading good samples from a good blockette
                //--------------------------------------------
                DecompressedData decompData = currentRecord.decompress();
                float[] temp = decompData.getAsFloat();
                for (int i =0; i< temp.length; i++) realsamples.add((float)temp[i]);
                lgs=temp[temp.length-1];
                prevRecord = currentRecord;
            }

            if (!realsamples.isEmpty()){
                //--------------------------------------------
                // ArrayList<Float> ----> float[]
                //--------------------------------------------
                float[] data = new float[realsamples.size()];
                for (int i = 0; i < realsamples.size(); i++) {
                    data[i] = realsamples.get(i);
                }

                return data;
            } else return null;

        } catch (Exception ex) {

            return null;
        }

    }