Mismatch between image size and meta-data field of view

[pless]

Description

Field of view in metadata is not correct and is causing issues with subestting.

This is the original issue as reported by @pless :

In the Stereo image data, the json file has two fields that may be inconsistent with each other:

in:

"sensor_fixed_metadata": {...
      "field of view at 2m in X- Y- direction [m]": "[1.857 1.246]"... }

"sensor_variable_metadata": {...
     "height Left Image [pixel]": "2472",
      "width Left Image [pixel]": "3296"}

These are inconsistent because the field of view is a 3 to 2 aspect ratio and the image size is a 4 to 3 aspect ration.  While it is possible that the camera has pixels that are not square, this is unusual.  Perhaps the calculation for the field of view was based on the x-axis and the conversion to the y-axis field of view just assumed a 3-2 ratio?

Context

This field of view estimate is used on the creation of the geoTIFFs and the overall field summaries. The version created for the June 2 event required "magic numbers" to fix the scale and this may have been a partial cause of that.

Completion Criteria

• [ ] Check that field of view meta-data in sensor metadata is correct
• [ ] Update extractors to query values from metadata in Clowder dynamically

[smarshall-bmr]

I'm guessing that the field of view is incorrect. Allied Vision does indeed list the sensor as a 4/3 with square 5.5 µm x 5.5 µm pixels. https://www.alliedvision.com/en/products/cameras/detail/Prosilica%20GT/3300.html

[ghost]

@smarshall-bmr - can you provide us with the correct field of view? Priority is stereo-top

[max-zilla]

you can see examples of the misalignments that result from lack of FOV consideration in #265

[smarshall-bmr]

I'm dropping hyperspectral for a moment to do this, check back soon.

[smarshall-bmr]

At 2m the field of view of each camera is 101.5x74.9cm or 28.4 degrees on the x and 21.2 degrees on the y axis if my trig serves me correctly.

[smarshall-bmr]

@pless @max-zilla @craig-willis @rachelshekar @dlebauer

Just to prevent any panic about how different these values are from the metadata reported ones I wanted to say that I roughly checked the FOV before programming StereoVIS scans. The current scans are based on a 100cmx75cm estimated FOV so scans take images every 50cm on the y axis and the gantry moves 105cm on the x between scans with a total of 3 scans per range. The StereoVIS cameras are 30cm apart so a single measurement actually captures around 130cm on the x-axis at 2m. This means each range should have around 270cm of 2x coverage and around 340cm total coverage per scan on 350cm ranges.

I'm certainly welcome to changes to this scanning procedure if anyone has suggestions.

[max-zilla]

@smarshall-bmr thanks. @ZongyangLi I am going to run some new results with this FOV:

"field of view at 2m in X- Y- direction [m]": 

OLD:
1.857 x 1.246
NEW:
1.015 x 0.749

[max-zilla]

@smarshall-bmr @ZongyangLi @yanliu-chn I added those new FOV to the get_fov() method of bin_to_geotiff():

def get_fov(metadata, camHeight, shape):
    try:
        cam_meta = metadata['lemnatec_measurement_metadata']['sensor_fixed_metadata']
        fov = cam_meta["field of view at 2m in x- y- direction [m]"]
    except KeyError as err:
        fail('Metadata file missing key: ' + err.args[0])

    try:
        fov_list = fov.replace("[","").replace("]","").split()
        fov_x = float(fov_list[0])
        fov_y = float(fov_list[1])

        # TODO: Find a better solution once metadata files are fixed
        # TODO: These values from https://github.com/terraref/computing-pipeline/issues/126#issuecomment-292027575
        fov_x = 1.015
        fov_y = 0.749

        # test by Baker
        gantry_meta = metadata['lemnatec_measurement_metadata']['gantry_system_variable_metadata']
        gantry_z = gantry_meta["position z [m]"]
        fov_offset = (float(gantry_z) - 2) * FOV_MAGIC_NUMBER
        fov_y = fov_y*(FOV_IN_2_METER_PARAM + fov_offset)
        fov_x = (fov_y)/shape[1]*shape[0]

        # given fov is at 2m, so need to convert for actual height
        #fov_x = (camHeight * (fov_x/2))/2
        #fov_y = (camHeight * (fov_y/2))/2

...but that is not sufficient:

Left  - stereoTop - 2017-02-09__12-52-49-893 (right TIF)
Right - stereoTop - 2017-02-09__12-52-48-383 (right TIF)

...you can see the little half-crescent shadow that does not align. There are also alignment issues between right/left TIF: Here, I've turned on the left halves of the stereo image (which appear shifted up from the right half in the first image) - see the crescent shadow again, offset too far up.

At this point I want to turn fixing this problem over to a more knowledgable party. I think it would be best to use these two datasets as a small sample and get the demosaic bin_to_geotiff.py script working properly with these before we proceed.

Links to datasets in Clowder: https://terraref.ncsa.illinois.edu/clowder/datasets/58de9a844f0c430e2c3fb56f https://terraref.ncsa.illinois.edu/clowder/datasets/58de9a934f0c430e2c3fb823

1) download both datasets 2) use the .bin files as basis for testing 3) use QGIS (free) for mapping and checking georegistration

Who is best suited to doing this? @yanliu-chn do you know of someone in CyberGIS with time to troubleshoot this? @ZongyangLi how comfortable are you with this aspect of the processing? In particular I wonder about the FOV offset portion of the code block I pasted above.

It may not be as simple as fixing the FOV parameters, but this is a critical piece for the rest of the pipeline to get fixed.

[smarshall-bmr]

@max-zilla Do you think we could use the initial georegistration as a starting point for a feature matching algorithm to do final alignment?

[ZongyangLi]

@max-zilla I am using the following code to estimate the FOV now.

def get_fov(metadata, camHeight, shape): fov_x = 1.015 fov_y = 0.749 HEIGHT_MAGIC_NUMBER = 1.64 PREDICT_MAGIC_SLOPE = 0.574 predict_plant_height = PREDICT_MAGIC_SLOPE * camHeight camH_fix = camHeight + HEIGHT_MAGIC_NUMBER - predict_plant_height fix_fov_x = fov_x*(camH_fix/2) fix_fov_y = fov_y*(camH_fix/2) return (fix_fov_x, fix_fov_y)

This code base on two assumption:

1. cameraHeight(m)/2(m) = fov_camera/fov_in_2m
2. ground may not visible because of the canopy cover

But I still need to add these two 'HEIGHT_MAGIC_NUMBER' and 'PREDICT_MAGIC_SLOPE' to my formula to let the stitching things looks reasonable for all the stereoTop dataset.

[ZongyangLi]

Stitching map looks fine to me both for ground visible and ground non-visible.

[max-zilla]

@ZongyangLi can you list which datasets/images you're using to show those examples?

[ZongyangLi]

@max-zilla They are datasets among '2016-05-0712-59-33-639' and '2016-05-0713-01-03-699', '2016-06-3009-23-46-258' and '2016-06-3009-24-02-948'.

[max-zilla]

@ZongyangLi can you try a few spot-checkes with 2017 recent data? I ran the pipeline on some data from 02-10:

[ZongyangLi]

@max-zilla Here is the stitching result in those images.

If this works for you, I could update the code to the repo.

[dlebauer]

@ZongyangLi could you please write a test that passes if and only if the images are properly aligned?

[dlebauer]

@pless could you please review this issue and comment on what a test for proper alignment would look like? What is good enough? What are appropriate input files and output conditions?

[solmazhajmohammadi]

@max-zilla @ZongyangLi If the mismatch happened in the recent dataset, It might be caused by the wheel change back in November 22, 2016

[max-zilla]

The mismatches I've seen have all been from Feb 2017.

[ZongyangLi]

@max-zilla @dlebauer Maybe I didn't describe clear enough here: https://github.com/terraref/computing-pipeline/issues/126#issuecomment-293368555 I think you may not get a properly aligned map without update the code I mentioned. I have updated my code here: https://github.com/terraref/computing-pipeline/tree/demosaic_extractor/scripts/stereoImager

With the code 'full_day_to_tiles.py', I can get a well aligned map for almost all the stereoTop images. Or if that's not works for you, we may discuss with others to find a solution. By using 'full_day_to_tiles.py', you could pass 'python full_day_to_tiles.py -i input_parent_dir -o output_dir -d child_dir' to create the out put in 'output_dir', here 'input_parent_dir' means a directory as '/ua-mac/raw_data/stereoTop/', and 'child_dir' means whatever directory in the parent dir you want to process, like '2017-02-10'.

[max-zilla]

@ZongyangLi your comment made sense, I had implemented the code you posted on the extractor and ran it to get the latest images I posted. But perhaps there is other updated code I was missing.

I will check your updated code... note that we are not using the directory you linked anymore! I created a new repo for this code several months ago: https://github.com/terraref/extractors-stereo-rgb/tree/master/demosaic This is the repo that should be used and updated in the future. I'll make sure your latest code from here: https://github.com/terraref/computing-pipeline/tree/demosaic_extractor/scripts/stereoImager ...is included. From now on you can make pull requests against the stereo-rgb repo.

However I want to be clear that I am not the right person to decide if this is being done properly :) I will generate new results today but between @pless @dlebauer @yanliu-chn @robkooper we should make sure the right GIS folks are examining our results before saying it's good.

Thanks for pushing on this.

[max-zilla]

@ZongyangLi @pless these results are indeed looking better.

4:25 only:

4:25 + 4:24: ^^ you can still see a tiny amount of duplication in the middle of the tire tread, but this looks very close. I'll let this run and we can use this as an evaluation basis.

The script is currently running on 02-11: https://terraref.ncsa.illinois.edu/clowder/collection/58de9ce74f0c430e2c402d69

[ZongyangLi]

@max-zilla The magic number I tested comes from left side image, It will be closer when you use left side image. I didn't try right image yet because I am not sure whether the method I used is a good way to solve the problem.

[ghost]

David suggests that offset be annotated so users know that we are aware of it.

[pless]

I think this has gotten to a good point. There remains a magic number in the process which I don't understand (but must have to do with some coordinate system scaling), but Zongyang has derived a way to automatically reset this magic number as the gantry height changes.

@rachelshekar, I think that the correct annotation to this should say something like:

"This stitched image is computed based on an assumption that the scene is planar. There are likely to be be small offsets near the boundary of two images anytime there are plants at the boundary (because those plants are higher than the ground plane), or where the dirt is slightly higher or lower than average."

[max-zilla]

Excellent. Will get this started up. I think for clarity I am going to create a new Level_1 output going forward:

/ua-mac/Level_1/stereoTop_fullField 

...that will have the full field result. This will operate on another new Level_1 directory that will include the raw -> geoTiff temporary files:

/ua-mac/Level_1/stereoTop_geotiff

So pipeline will be: 1) datasets come in every few seconds, are converted to TIF and JPG in the 2nd directory 2) full field stitching extractor will be executed using Rule Checker to look for an entire day's worth of geotiffs once they are done in the 2nd directory 3) full field output is written to the 1st directory 4) later, eventually, we could clean out the _geotiffs directory and just keep the full field.

[jdmaloney]

@max-zilla I'd say create a directory in /gpfs/largeblockFS/scratch/terraref/stereoTop_geotiff (on nodes their should be a symlink from /scratch to /gpfs/largeblockFS/scratch, so you can use that shortened path if you want) This way the ephemeral data doesn't live within the production data structure.

[max-zilla]

@ZongyangLi in your full_day_to_tiles.py, I have a question:

# provide an output folder
        print("tiflist: %s" % tifflist)
        temp_out_dir = tifflist.replace(os.path.basename(tifflist), "")
        #test folder
        #temp_out_dir = '/Users/nijiang/Desktop/pythonTest/clowder3/'
        if not os.path.exists(temp_out_dir):
            os.makedirs(temp_out_dir)
        print("base dir: %s" % temp_out_dir)
        upper_dir = os.path.abspath(os.path.join(os.path.dirname(temp_out_dir), os.path.pardir))

        # Create VRT from every GeoTIFF
        print "Starting VRT creation..."
        full_day_to_tiles.createVrt(temp_out_dir,tifflist)
        print "Completed VRT creation..."

!!!!!!!!!!!

        # Generate tiles from VRT
        print "Starting map tile creation..."
        full_day_to_tiles.createMapTiles(temp_out_dir,multiprocessing.cpu_count())
        print "Completed map tile creation..."

        # Generate google map html template
        print "Starting google map html creation..."
        full_day_to_tiles.generate_googlemaps(temp_out_dir)
        print "Completed google map html creation..."

        print("Uploading output html to dataset")
        #html_out = os.path.join(temp_out_dir, 'opengooglemaps.html')

This looks like it stitches the GeoTIFFs into a large VRT, then chops the VRT into google maps-ready tiles.

For our extractor, we basically just want to stop at the !!!!!!!! for now and write that VRT to the permanent full-field output, right? The google tiles part we can do separately if we want.

[ZongyangLi]

@max-zilla I am not quite sure if we can stop before 'Generate tiles from VRT'.

The .vrt is a virtual TIFF, which is a way of storing all of the information for the overlapping TIFFs without actually creating the full stitched TIFF. Generating that full TIFF would result in an absurdly large file, so having this "virtual TIF" is a really nice alternative, and one that the map tiler handles well. The final piece is that we then generate the map tiles from the VRT.

So if you have some other way of showing stitching map using the information in VRT or it would be better to do that separately when we want, that will be fine. Let me know if this dose not make sense.

[max-zilla]

@ZongyangLi the goal of this extractor is to generate an absurdly large TIFF, I think... can you load a VRT into GIS? Can you view it like a TIF file? The new plan is to generate the plot-level clips of the large VRT on-demand, not generating them all every time automatically.

If users are happy with a VRT, my thinking was to load the VRT into Clowder and let people use that. Is the VRT just a metadata file of sorts that points to all the component TIFF images?

@jterstriep @yanliu-chn @dlebauer @robkooper do you have experience with VRTs for the use we want?

[ZongyangLi]

@max-zilla To my understanding, this VRT just contains tif path and Geo-information to the tif file. to view the large stitching map, you have to create tile files.

For our new plan, to generate the plot-level clips, we may first use the metadata to collect all the 'in the plot' tiffs, then use these classified tiffs to create VRT and plot-level stitching map.

[yanliu-chn]

Yes, we handle a lot of VRTs.

VRT is open source geospatial software community "standard". Esri has a similar format called mosaic dataset. As Zongyang said, it is just an XML file that contains geo metadata such as projection, extent, resolution, plus a list of file entries, each including attributes of a file included in VRT, e.g., data type, spatial metadata, file path, and relative position in the VRT extent. It is a very convenient way to handle large or a large number of files.

It is important to note that: if you publish VRT as a data product, you have to make sure all the consuming clients of the VRT know how to interpret the file path. If you put the absolute path on ROGER as file path, that path has to be accessible by all the clients, be it a VM mounting ROGER storage as an NFS mount point or a desktop software. So I highly recommend that you use relative path. There are two ways to use relative path. One is to directly put a path without the beginning "/"; the other is to use the "relativeToVRT" option in XML file entry. Either is fine.

Of course, then, you have to make sure that that relative path never changes; otherwise, VRT needs to be updated accordingly.

Most of GIS software understands VRT, so as ArcGIS. You can treat VRT the same as GeoTIFF in that regard.

If files in a VRT have overlapping spatial extent, it is an unknown behavior how the overlapping is handled. VRT doesn't handle it. GDAL does by picking the value of the overlapped area from one (either the first or the last) file. For areas that no file has content on, most of GIS software returns nodata value.

[yanliu-chn]

One example of a VRT file. This is the VRT for the national elevation dataset at 10m resolution. We use absolute path in it.

<VRTDataset rasterXSize="637212" rasterYSize="280812">
  <SRS>GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]]</SRS>
  <GeoTransform> -1.2500055555560000e+02,  9.2592592593335375e-05,  0.0000000000000000e+00,  5.0000555555555515e+01,  0.0000000000000000e+00, -9.2592592593335375e-05</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-3.402823466385289e+38</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/gpfs_scratch/usgs/ned10m/n25w081.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="10812" RasterYSize="10812" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="10812" ySize="10812" />
      <DstRect xOff="475200" yOff="270000" xSize="10812" ySize="10812" />
      <NODATA>-3.402823466385289e+38</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/gpfs_scratch/usgs/ned10m/n25w082.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="10812" RasterYSize="10812" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="10812" ySize="10812" />
      <DstRect xOff="464400" yOff="270000" xSize="10812" ySize="10812" />
      <NODATA>-3.402823466385289e+38</NODATA>
    </ComplexSource>
... ...

[max-zilla]

@czender @ZongyangLi @yanliu-chn @dlebauer @jterstriep I have generated a VRT for 2017-04-27...

/sites/ua-mac/Level_1/fullfield/2017-04-27/stereoTop_fullfield.VRT

the VRT file is only 13 MB.

I'm manually running

gdal_translate stereoTop_fullfield.VRT stereoTop_fullfield.tif

in that same directory on the VM to see what kind of size/look the output TIF has. but it's gonna be massive and take a while - after 45 mins this morning, it's at like 2%.

[yanliu-chn]

I'm not sure how useful the big geotiff is. Even ppl are patient to download it, desktop software would have hard time to display it.

[max-zilla]

@yanliu-chn is there a simple arg to gdal_translate to downsample that might speed up generation and still give a useful "human-readable" idea of the data? The VRT still requires you to download 8,000 (!) geoTIFFs in this case or have visibility to them on Roger.

[yanliu-chn]

@max-zilla yes. use the "-outsize" option to give it a percentage or num of pixels for the output geotiff. And select a good resampling method with "-r". I don't know which resampling method is better for such images. @czender or @ZongyangLi can decide. Nearest neighbor is the fastest, but we rarely use it because it is too simple.

[max-zilla]

thanks. it's gone to 10% in the last half hour so I'm gonna wait and see if it finishes by the end of day today.

[ZongyangLi]

This mismatch issue in stereo camera can be closed. Further discussion #265