search

davidfrantz / force

Framework for Operational Radiometric Correction for Environmental monitoring
GNU General Public License v3.0
172 stars 50 forks source link

Coordinate handling in force-level1-csd #136

Closed davidfrantz closed 2 years ago

davidfrantz commented 2 years ago

I just stumbled upon some issues in force-level1-csd:

When naively using a coordinate, no images are found:

force-level1-csd -n -s S2A -c 0,70 . . queue.txt "7/50"

Searching for footprints / tiles intersecting with input geometry...

Querying the metadata catalogue for Sentinel2 data
Sensor(s): S2A
Tile(s): PRFID,
Daterange: 19700101 to 20211028
Cloud cover minimum: 0%, maximum: 70%

There were no Sentinel2 Level 1 scenes found matching the search criteria.

Done.

Probably my bad as this is not supported? But an error message would be nice. When using a bounding box, it works as expected for Sentinel-2:

force-level1-csd -n -s S2A -c 0,70 . . queue.txt "7/50,7/51,8/51,8/50,7/50"

Searching for footprints / tiles intersecting with input geometry...

Querying the metadata catalogue for Sentinel2 data
Sensor(s): S2A
Tile(s): T31UGR,T31UGS,T32ULA,T32ULB,T32UMA,T32UMB
Daterange: 19700101 to 20211028
Cloud cover minimum: 0%, maximum: 70%

1860 Sentinel2 Level 1 scenes matching criteria found
837.81GB data volume found.

Done.

When switching to Landsat, it fails though:

force-level1-csd -n -s LC08 -c 0,70 . . queue.txt "7/50,7/51,8/51,8/50,7/50"

Searching for footprints / tiles intersecting with input geometry...

Querying the metadata catalogue for Landsat data
Sensor(s): LC08
Tier(s): T1
Tile(s): "196024","196025","197024","197025","195025"
Daterange: 19700101 to 20211028
Cloud cover minimum: 0%, maximum: 70%

There were no Landsat Level 1 scenes found matching the search criteria.

Done.

The Path/Rows are correct though. When using one of these directly, it works again:

force-level1-csd -n -s LC08 -c 0,70 . . queue.txt 196025

Querying the metadata catalogue for Landsat data
Sensor(s): LC08
Tier(s): T1
Tile(s): 196025
Daterange: 19700101 to 20211028
Cloud cover minimum: 0%, maximum: 70%

84 Landsat Level 1 scenes matching criteria found
83.25GB data volume found.

Done.

Another issue: I believe that the out-of-bounds checking of the coordinates is done on the wrong axis (using a longitude of 100):

force-level1-csd -n -s LC08 -c 0,70 . . queue.txt "100/50,7/51,8/51,8/50,7/50"

force-level1-csd - FORCE cloud storage downloader for Landsat and Sentinel-2
https://force-eo.readthedocs.io/en/latest/components/lower-level/level1/level1-csd.html

Usage: force-level1-csd [-c min,max] [-d starttime,endtime] [-n] [-k] [-s sensor]
                        [-t tier] [-u]
                        metadata-dir level-1-datapool queue aoi

Error: Latitude out of range
       Coordinate: 100/50 - 100 is not in range -90 to 90
       This error may also mean that you tried to use a vector file as AOI but provided an incorrect path
ernstste commented 2 years ago

Hi David,

  1. Coordinate input is exclusively for polygons, as stated in the documentation. I agree that a check for too few coordinates could make sense though. Or even switching to creating a point layer if there is only one coordinate provided.
  2. Looks like the path/rows got casted in quotes in the Landsat query. I'm not sure why this is the case. Just tested it here and it happened on R2D2 but not on BB8 on 3.7.2.
  3. Looks like you're right. Just glanced a the code and the variable assignment was mixed up.

Stefan

davidfrantz commented 2 years ago

Thanks Stefan,

  1. looks like it is OS-related. BB8 still runs on Ubuntu 18.04. R2D2 and the Docker (that I used) on 20.04.

Cheers, David

ernstste commented 2 years ago

Interesting. I'm not sure why it would only happen for S2 though, the data types of the corresponding field are the same for LS and S2. Here is the command, returning the results casted for LS: ogr2ogr -f CSV /vsistdout/ -select PRFID 'WFS:http://ows.geo.hu-berlin.de/cgi-bin/qgis_mapserv.fcgi?MAP=/owsprojects/grids.qgs&SERVICE=WFS&REQUEST=GetFeature&typename=sentinel2&Filter=%3Cogc:Filter%3E%3Cogc:Intersects%3E%3Cogc:PropertyName%3Eshape%3C/ogc:PropertyName%3E%3Cgml:Polygon%20srsName=%22EPSG:4326%22%3E%3Cgml:outerBoundaryIs%3E%3Cgml:LinearRing%3E%3Cgml:coordinates%3E7,50%207,51%208,51%208,50%207,50%3C/gml:coordinates%3E%3C/gml:LinearRing%3E%3C/gml:outerBoundaryIs%3E%3C/gml:Polygon%3E%3C/ogc:Intersects%3E%3C/ogc:Filter%3E'

ernstste commented 2 years ago

I opened a PR which fixes the range check of the coordinates and adds the possibility of using single coordinates representing a point as AOI. Haven't found the time to look into the issue with the quoted results being returned by ogr2ogr though and I don't really have an idea what causes it at the moment.

davidfrantz commented 2 years ago

OK Thanks, I will look into the PR.

The ogr command above is for S2 actually. Can you post the one for Landsat, too?

davidfrantz commented 2 years ago

Got it ogr2ogr -f CSV /vsistdout/ -select PRFID 'WFS:http://ows.geo.hu-berlin.de/cgi-bin/qgis_mapserv.fcgi?MAP=/owsprojects/grids.qgs&SERVICE=WFS&REQUEST=GetFeature&typename=landsat&Filter=%3Cogc:Filter%3E%3Cogc:Intersects%3E%3Cogc:PropertyName%3Eshape%3C/ogc:PropertyName%3E%3Cgml:Polygon%20srsName=%22EPSG:4326%22%3E%3Cgml:outerBoundaryIs%3E%3Cgml:LinearRing%3E%3Cgml:coordinates%3E7,50%207,51%208,51%208,50%207,50%3C/gml:coordinates%3E%3C/gml:LinearRing%3E%3C/gml:outerBoundaryIs%3E%3C/gml:Polygon%3E%3C/ogc:Intersects%3E%3C/ogc:Filter%3E'

davidfrantz commented 2 years ago

This would fix it:

ogr2ogr -f CSV /vsistdout/ -select PRFID 'WFS:http://ows.geo.hu-berlin.de/cgi-bin/qgis_mapserv.fcgi?MAP=/owsprojects/grids.qgs&SERVICE=WFS&REQUEST=GetFeature&typename=landsat&Filter=%3Cogc:Filter%3E%3Cogc:Intersects%3E%3Cogc:PropertyName%3Eshape%3C/ogc:PropertyName%3E%3Cgml:Polygon%20srsName=%22EPSG:4326%22%3E%3Cgml:outerBoundaryIs%3E%3Cgml:LinearRing%3E%3Cgml:coordinates%3E7,50%207,51%208,51%208,50%207,50%3C/gml:coordinates%3E%3C/gml:LinearRing%3E%3C/gml:outerBoundaryIs%3E%3C/gml:Polygon%3E%3C/ogc:Intersects%3E%3C/ogc:Filter%3E' | sed 's/"//g'

ernstste commented 2 years ago

Hmmm I would prefer knowing what's going on there to not run into more issues further down. But I guess it works as a quick fix.

ernstste commented 2 years ago

I pushed the suggested fix. It's still a mystery to me though. The same does not happen when using geopackage in the ogr2ogr call.

davidfrantz commented 2 years ago

Thanks Stefan, I just pulled your changes! Cheers, David

punyidea commented 2 years ago

I think this might need to be looked into some more. I tried running the program with this simple bounding box for Munich and got tiles with flipped coordinates.

force-level1-csd -c 0,10 -k -s lc08 metadata test_data/level1 test_data/l1_pool.txt city_bb.gpkg (sorry, can't upload geopackage)

Bounding box command also has mistake. force-level1-csd -n -s lc08 metadata test_data/level1 test_data/l1_pool.txt 48.139/11.547,48.150/11.547,48.150/11.574,48.139/11.574,48.139/11.547

punyidea commented 2 years ago

Hello everyone, please disregard the above. I realize that simply replacing the bash file with update was not sufficient- after updating accordingly I did not experience the coordinate bug. I tested the coordinate version and it worked, after also seeing coordinates were incorrectly put.

However, I still have the same issue with the Geopackage file - transposing coordinates.

ernstste commented 2 years ago

Hi punyidea,

The update did not change the coordinate handling itself, but only the check that verifies if the coordinates are in the correct range. So the results of your query would not change.

Could you clarify what you mean when you say 'transposing coordinates'? What is the ouptut you expect and what is the output you receive?

punyidea commented 2 years ago

Hi enstste,

Thank you for your quick reply.I'll describe my workflow so that it can be replicated.

1) I created and saved a simple polygon bounding box, of e.g. the city of Munich using qgis (see German Language screenshot of export). grafik

2) I used the bounding box as an AOI. 3) The tiles returned by the tile query correspond to coordinates with flipped latitude and longitude, in this case over the Arabian peninsula and Horn of Africa. grafik

ernstste commented 2 years ago

Thanks! I was able to reproduce the issue. It seems like the Landsat path/row features that are downloaded from the server are incorrect for some reason. The piece of code that is responsible for pulling them has not changed. Will have to see what's going on here.

davidfrantz commented 2 years ago

To add here: we just ran into this problem, too. Both Landsat and Sentinel-2 were downloaded for Western Sahara instead of Zambia when using a gpkg.

davidfrantz commented 2 years ago

Hi @ernstste, could you please have another look at this? Many thanks, David

ernstste commented 2 years ago

Hi David, I haven't forgotten about the issue. Sorry, but this is simply beyond my capacity at the moment. I do believe we should open the discussion about shipping force with the tiling grids for s2 and landsat again to avoid running into weird issues like this one or the last one. We're not really gaining anything from the added complexity that comes with querying the data from a server every time.

punyidea commented 2 years ago

I would also like to comment that I found it simpler in the end to download Landsat data using the EarthExplorer interface, which had all of the features that I personally needed from the downloader. Cloud cover and time frame are included there as well.

davidfrantz commented 2 years ago

I do believe we should open the discussion about shipping force with the tiling grids for s2 and landsat again to avoid running into weird issues like this one or the last one. We're not really gaining anything from the added complexity that comes with querying the data from a server every time.

I believe it would be good to not have the grids in the software repository as the vector files are fairly large. But what about uploading them to a long-term data repository like Zenodo? Would that help?

davidfrantz commented 2 years ago

I would also like to comment that I found it simpler in the end to download Landsat data using the EarthExplorer interface, which had all of the features that I personally needed from the downloader. Cloud cover and time frame are included there as well.

Has this become fast and stable enough for larger areas and time series?

ernstste commented 2 years ago

I believe it would be good to not have the grids in the software repository as the vector files are fairly large. But what about uploading them to a long-term data repository like Zenodo? Would that help?

Good idea, that could be a viable alternative: download the files if they are not already present in the same directory as the metadata catalogue. I should have much more time from next month and will look into possible solutions

punyidea commented 2 years ago

I would also like to comment that I found it simpler in the end to download Landsat data using the EarthExplorer interface, which had all of the features that I personally needed from the downloader. Cloud cover and time frame are included there as well.

Has this become fast and stable enough for larger areas and time series?

I used the Bulk Downloader interface that they recommend, and downloaded timeseries level 1 and 2 data for two separate regions containing cities without much issue. It has a downside that (at least I have not tried otherwise) I navigated their website to download and select which files to use, and so it's not as simple as a simple command run from the terminal.

Overall my experience here is a bit limited but I would certainly recommend trying it out - the learning curve here is not steep. https://earthexplorer.usgs.gov/

davidfrantz commented 2 years ago

I would also like to comment that I found it simpler in the end to download Landsat data using the EarthExplorer interface, which had all of the features that I personally needed from the downloader. Cloud cover and time frame are included there as well.

Has this become fast and stable enough for larger areas and time series?

I used the Bulk Downloader interface that they recommend, and downloaded timeseries level 1 and 2 data for two separate regions containing cities without much issue. It has a downside that (at least I have not tried otherwise) I navigated their website to download and select which files to use, and so it's not as simple as a simple command run from the terminal.

Overall my experience here is a bit limited but I would certainly recommend trying it out - the learning curve here is not steep. https://earthexplorer.usgs.gov/

I have stopped using this tool several years ago because it was slow and unstable - especially when downloading a lot of data. But I am happy to hear that this has become a viable alternative again.

davidfrantz commented 2 years ago

Hi @ernstste,

this issue here is actually solved, right?

Cheers, David

ernstste commented 2 years ago

It is. Sorry I thought I had set the PR to close this issue.