Queries regarding force-tile-extent and associated issues

[paulvpop]

I went to the documentation for force-tile-extent (it seems to be old. More on this below), and for this a datacube-definition.prj was needed (as shown here). So, I made one by writing the following and saving it using unicode style file:

PROJCS["WGS 84 / India NSF LCC",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",24],PARAMETER["central_meridian",80],PARAMETER["standard_parallel_1",12.472955],PARAMETER["standard_parallel_2",35.1728044444444],PARAMETER["false_easting",4000000],PARAMETER["false_northing",4000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","7755"]]
93.9574520
29.3477290
49.2859027
-14.0282153
30000.000000
3000.0000000

I have faced a few issues. I will explain them in detail:

The projection I chose was EPSG:7755.

When I run the command davidfrantz/force gdalsrsinfo -e EPSG:7755 I get the details correctly and when I run davidfrantz/force gdalsrsinfo -V -e EPSG:7755 , I get the "Validate Succeeds" message with the details.

However, when I Copy Text the wkt from https://epsg.io/7755 to get it in one line, and run the davidfrantz/force gdalsrsinfo -V 'PROJCS["WGS 84 / India NSF LCC",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",24],PARAMETER["central_meridian",80],PARAMETER["standard_parallel_1",12.472955],PARAMETER["standard_parallel_2",35.1728044444444],PARAMETER["false_easting",4000000],PARAMETER["false_northing",4000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","7755"]]' I get the following error:

"ERROR 1: missing [
ERROR 1: ERROR - failed to load SRS definition from PROJCS[WGS 84 / India NSF LCC,GEOGCS[WGS 84,DATUM[WGS_1984,SPHEROID[WGS 84,6378137,298.257223563]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433,AUTHORITY[EPSG,9122]],AUTHORITY[EPSG,4326]],PROJECTION[Lambert_Conformal_Conic_2SP],PARAMETER[latitude_of_origin,24],PARAMETER[central_meridian,80],PARAMETER[standard_parallel_1,12.472955],PARAMETER[standard_parallel_2,35.1728044444444],PARAMETER[false_easting,4000000],PARAMETER[false_northing,4000000],UNIT[metre,1],AXIS[Easting,EAST],AXIS[Northing,NORTH],AUTHORITY[EPSG,7755]]"

But that doesn't make sense. I validated the wkt in QGIS too. Nothing's missing. Do you know why this happens?

I had created a bounding box for my study area to extract the north-west point's geographic longitude and latitude (in WGS84 format - 93.9574520
29.3477290 respectively) and the projected longitude (X) and latitude (Y) (in EPSG:7755 format - 49.2859027
-14.0282153 respectively). Is this correct? I am asking because I am a little confused by the custom projection you had used in this example. 2456026.250000,
4574919.500000 are longitude and latitude of the projected coordinate system?

3. The documentation for force-tile-extent seems old, but the documentation in program itself is updated. I get the following when I run davidfrantz/force force-tile-extent -h:

"Usage: force-tile-extent [-hvi] [-d datacube-dir] [-a allow-list] [-j jobs] input-vector

optional:
-h = show this help
-v = show version
-i = show program's purpose

-d = directory of existing datacube
defaults to current directory
'datacube-definition.prj' needs to exist in there
-a = tile allow-list to restrict the processing extent
The extent will be computed based on the input AOI
Defaults to './allow-list.txt'
-j = number of jobs, defaults to 100%

Positional arguments:

• input-vector: a polygon vector file for the area of interest (AOI)

---

see https://force-eo.readthedocs.io/en/latest/components/auxilliary/tile-extent.html"

So, I modified the code to force-tile-extent -d /home/docker/force/force/datacube-dir -a /home/docker/force/force/allow-list.txt -j 100% /home/docker/force/force/Siang_SA_Prj/Siang_SA_Prj.shp

and got the following ouput:

Suggested Processing extent parameters:
X_TILE_RANGE = 178 183
Y_TILE_RANGE = -156 -149
FILE_TILE = /home/docker/force/force/allow-list.txt

Maybe, you can give this as an example for your documentation page.

4. I am getting negative tile numbers. Here's the file output of the above command:

48
X0178_Y-149
X0178_Y-150
X0178_Y-151
X0178_Y-152
X0178_Y-153
X0178_Y-154
X0178_Y-155
X0178_Y-156
X0179_Y-149
X0179_Y-150
X0179_Y-151
X0179_Y-152
X0179_Y-153
X0179_Y-154
X0179_Y-155
X0179_Y-156
X0180_Y-149
X0180_Y-150
X0180_Y-151
X0180_Y-152
X0180_Y-153
X0180_Y-154
X0180_Y-155
X0180_Y-156
X0181_Y-149
X0181_Y-150
X0181_Y-151
X0181_Y-152
X0181_Y-153
X0181_Y-154
X0181_Y-155
X0181_Y-156
X0182_Y-149
X0182_Y-150
X0182_Y-151
X0182_Y-152
X0182_Y-153
X0182_Y-154
X0182_Y-155
X0182_Y-156
X0183_Y-149
X0183_Y-150
X0183_Y-151
X0183_Y-152
X0183_Y-153
X0183_Y-154
X0183_Y-155
X0183_Y-156

So, I went through the documentation and it has already been addressed in the "Frequently Asked Questions"

"The tile IDs of the processed Level 2 data have negative numbers, e.g. X-100_Y0100.

Make sure that the origin of the target grid (ORIGIN LAT / ORIGIN LON) is in the North-West of your study area. Potentially, you have accidentally swapped latitude and longitude. Note that a geographic location in the North-West is not necessarily North-West in the output coordinate system, too (for an example see Fig. 14). Although not recommended, higher-level FORCE functions should be able to digest negative tile numbers (note that we did not test this exhaustively)."

I couldn't find Fig. 14. Could you share the link with me? I am trying to figure out why I am getting negative tile numbers even though I selected the north-west points of the bounding box (attaching it here)
Rectangle_bounding_Siang.zip
containing the study area in both geographic and projected coordinate systems, and made sure that the longitude and latitude are in the right place.

Any idea why I am still getting negative tile numbers?

5. Also, a related minor issue I found was there are things like VIII.M and VII.X mentioned in some of the documentation pages. For example "Any gridded data cube (containing a data cube definition file, see VII.M) can be given as 1st argument.". What do these roman numerals refer to?

Thanks in advance for all the help.

With kind regards,
Paul

Originally posted in #350 (reply in thread)

[davidfrantz]

Hi @paulvpop,

thanks for opening a new thread. There are multiple sub-questions here. I try to break it down.

1. yes, it seems that the documentation of force-tile-extent is outdated. I have this on my to-do list, though on a lower priority for the moment.

2. If you want to have an overview of different programs, best use force-info. This is not 100% finished, but still gives you a good overview of what is there, e.g. this one: force-cube-init: Initialize a datacube definition. This tool allows you to generate a datacube-definition.prj based on a Level 2 parameter file without actually running the preprocessing before.

3. this might be related to how parameters are passed through Docker. If I use plain Docker, it works on my end:
   ``docker run davidfrantz/force gdalsrsinfo -V 'PROJCS["WGS 84 / India NSF LCC",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",24],PARAMETER["central_meridian",80],PARAMETER["standard_parallel_1",12.472955],PARAMETER["standard_parallel_2",35.1728044444444],PARAMETER["false_easting",4000000],PARAMETER["false_northing",4000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","7755"]]'

Validate Succeeds

PROJ.4 : +proj=lcc +lat_0=24 +lon_0=80 +lat_1=12.472955 +lat_2=35.1728044444444 +x_0=4000000 +y_0=4000000 +datum=WGS84 +units=m +no_defs

OGC WKT2:2018 :
PROJCRS["WGS 84 / India NSF LCC",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["unnamed",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",24,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",80,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",12.472955,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",35.1728044444444,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",4000000,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",4000000,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",7755]]
``

4. regarding the negative tile numbers, please try force-cube-init first. Probably, the manual conversion between different coordinate systems was flawed.

5. the FAQ page should actually not be included anymore since it is way outdated. I just deleted it.

[paulvpop]

Thanks for the reply.

1. Got it.
2. Thanks. Will try it out.
3. Will trying debugging. I am using Podman instead of Docker.
4. Will try it out. Thanks.
5. Oh right. I hope all the documentation is the respective pages instead.