create grid and reproject

97 views
Skip to first unread message

moh713...@gmail.com

unread,
Oct 14, 2020, 12:49:11 PM10/14/20
to FORCE
Dear List,
as you can see the image attached, the study area(black dash line : it is an island) is covered by 4 sentinel-2 tiles. I have two inquires:
1- how can I create the ARD with less number of tiles (design the grid)
2- how can I edit the projection in the parameter file for the data cube creation?
 I hope this is very clear.
I use the latest version of FORCE on ubuntu 18.04
thanks in advance

image_study_area.png

Stefan

unread,
Oct 14, 2020, 2:35:01 PM10/14/20
to FORCE
Dear Mohammed,
this is well documented in the tutorials.
Stefan

moh713...@gmail.com

unread,
Oct 17, 2020, 5:43:50 PM10/17/20
to FORCE
Hello Stefan,
this is what really happens 
The data cube parameters :
# DATA CUBES
# ------------------------------------------------------------------------
DO_REPROJ = TRUE
DO_TILE = TRUE
FILE_TILE = /home/localuser/force/mis/allow-list.txt
TILE_SIZE = 30000
BLOCK_SIZE = 3000
RESOLUTION_LANDSAT = 30
RESOLUTION_SENTINEL2 = 10
ORIGIN_LON =  12     
ORIGIN_LAT =  53 
PROJECTION = PROJCRS["ProjWiz_Custom_Cylindrical_Equal_Area",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Cylindrical Equal Area (Spherical)",
            ID["EPSG",9834]],
        PARAMETER["Latitude of 1st standard parallel",6.3630421,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Longitude of natural origin",53.9222717,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

RESAMPLING = NN

and this is the log file after running
/home/localuser/force/level1/sentinel/1-2017/T39PZQ/S2A_MSIL1C_20170101T064252_N0204_R120_T39PZQ_20170101T064359.SAFE
: Transformation failed.
Computing tile origin in dst_srs failed. Starting datacube failed.
Starting datacube(s) failed.


so what is going wrong?
thanks in advance

david frantz

unread,
Oct 19, 2020, 7:16:36 AM10/19/20
to FORCE
Hi Mohammed, 

please check your settings for 
ORIGIN_LON =  12     
ORIGIN_LAT =  53 
I suspect that you have swapped X and Y.

Best,
David

moh713...@gmail.com

unread,
Oct 19, 2020, 11:02:55 AM10/19/20
to FORCE
Hi DAVID,
I have gone through the original code :

cube->origin_geo.x = pl2->orig_lon;
    cube->origin_geo.y = pl2->orig_lat;

    if ((warp_geo_to_any(cube->origin_geo.x,  cube->origin_geo.y,
                        &cube->origin_map.x, &cube->origin_map.y,
                         cube->proj)) == FAILURE){
      printf("Computing tile origin in dst_srs failed. ");
      free_multicube(multicube);
      return NULL;
    }
 trying to explore how I swapped X and Y but without getting progress to solve the problem.could you clarify it please?

moh713...@gmail.com

unread,
Oct 19, 2020, 11:09:52 AM10/19/20
to FORCE
hopefully with this image , the problem will be more clear
grid.png

david frantz

unread,
Oct 19, 2020, 11:11:49 AM10/19/20
to FORCE
Dear Mohammed, 
this has nothing to do with the code. Each projection is only valid within defined limits.
Your area of interest is at 12° Latitude and 53° Longitude.
In the parameter file, you defined 12°Longitude and 53° latitude (this is somewhere in Germany)
X = Longitude, Y = Latitude.

moh713...@gmail.com

unread,
Oct 19, 2020, 12:05:31 PM10/19/20
to FORCE
Dear David,
no things has been changed !!

david frantz

unread,
Oct 19, 2020, 12:40:35 PM10/19/20
to FORCE

1) Did you validate the projection?

2) if you cannot make your projection work, use one of the built-in projections, e.g. GLANCE7-AF

moh713...@gmail.com

unread,
Oct 19, 2020, 2:31:07 PM10/19/20
to FORCE
1- it is validated well
Validate Succeeds

PROJ.4 : +proj=cea +lat_ts=6.3630421 +lon_0=53.9222717 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

OGC WKT2:2018 :
PROJCRS["ProjWiz_Custom_Cylindrical_Equal_Area",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
2-   that already worked well but I want to use this projection to keep minimum distortion in area as  I plan to create land cover for this area .

hopefully this issue will be solved
thanks


moh713...@gmail.com

unread,
Oct 22, 2020, 9:09:38 AM10/22/20
to FORCE
Hi David,
After I validated the projection, I rewrote it in parameter file as below and it worked :
PROJECTION = PROJCRS["ProjWiz_Custom_Cylindrical_Equal_Area",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Lambert Cylindrical Equal Area (Spherical)",ID["EPSG",9834]],PARAMETER["Latitude of 1st standard parallel",6.3630421,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",53.9222717,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]

a new issue appeared when I checked the tiles as the most eastern part of the study area is not covered by the tiles as shown below :
force.png

I tried to check if the problem related to the projection choosing one from  the built-in projections but the problem is the same
thanks in advance.

david frantz

unread,
Oct 26, 2020, 10:53:49 AM10/26/20
to FORCE
Hi Mohammed,
I am not really sure what problem you are referring to.. Please elaborate more. From your image, I do not see any problem. Your study area is at the edge of a swath.. This means, the entire island will never be observed on the same day.
Best,
David

moh713...@gmail.com

unread,
Oct 26, 2020, 12:08:09 PM10/26/20
to FORCE
Hi David,
I used multiple images with different dates with a cloud cover less than 75%  covering all the study area(four tiles covering the whole  island as shown at the first message in this conversation).
the problem is that all tiles should  have to be produced to cover the whole study area as logically expected  and based on the input data, unfortunately the most right part of the study area has not been covered by the produced tiles as expected despite changing the projection, the cloud cover percentage of the processed images.

hopefully this is now very clear
Message has been deleted

Stefan

unread,
Oct 27, 2020, 6:01:39 AM10/27/20
to FORCE
Dear Mohammed,

see how the Sentinel-2 tiling grid relates to the Sentinel-2 orbits:

While FORCE provides a lot of easy to use tools, it does rely on its users having a solid base knowledge of remote sensing concepts and the particularities of the data used. In this free support group here we are trying to do our best to answer questions that are directly related to the functionality of FORCE. Providing support outside of this scope is unfortunately beyond our possibilites.

Stefan

moh713...@gmail.com

unread,
Oct 27, 2020, 6:16:35 AM10/27/20
to FORCE
Dear David, 
thank you for all the efforts you do for FORCE users.
I will do my best to see to which kind of problems this issue is related .
I will also try in another area to confirm from where this issue arise?
thanks again for  all your efforts  for the scientific community. 

moh713...@gmail.com

unread,
Oct 27, 2020, 6:27:34 AM10/27/20
to FORCE
Dear Stefan,
I think there is a misunderstanding to what I meant in my issue reporting.
I used in my process  sentinel-2 images that cover all the study area with characteristics that have already mentioned but after processing No tiles covering the right corner of the study area. I used all the multidate produced tiles to discover that issue.
I checked the cloud cover of the images to see if that is related to the parameters we set in the parameter file, I change the projection to the one in-built to see again but the issue is still .
hopefully this is very clear now

Stefan

unread,
Oct 27, 2020, 6:35:31 AM10/27/20
to FORCE
If you are 100% sure that you had L1 images from the orbit to the east, it's most likely that the L2 processing of these images failed due to high cloud cover / low registration targets. You can check the log files to see which scenes were processed successfully and which failed.

moh713...@gmail.com

unread,
Oct 27, 2020, 7:03:26 AM10/27/20
to FORCE
I really took into considerations all what mentioned and the log file did not report any issue regarding failing due to high cloud cover / low registration targets instead it reports this :
dc:  66.02%. wc:  66.14%. sc:   0.02%. cc:   0.36%. AOD: 0.1991. # of targets: 170/2.  1 product(s) written. Success! Processing time: 05 mins 57 secs

I will do the process again and see what is the issue. and report here.
thanks again

david frantz

unread,
Oct 27, 2020, 7:06:30 AM10/27/20
to FORCE
Dear Mohammed,
please indicate how many images you processed, and please post the original scene IDs.
Thanks,
David

moh713...@gmail.com

unread,
Oct 27, 2020, 7:24:51 AM10/27/20
to FORCE
Dear David,
Those are the original scene IDs : 
T39PYP
 T39PYQ
T39PZP
T39PZQ

Regarding the numbers of images I processed, I used four images covering the whole study area >
images with less 75 % cloud cover and  I tried with less than 5 % cloud cover when  I did not see the tiles covering the eastern part of the study area and the parameter file setting allow producing tiles to 75 % cloud cover.
when failed because the coregistration, I disable the coregistration in the parameter file but the problem is still.

Thanks

Stefan

unread,
Oct 27, 2020, 8:06:10 AM10/27/20
to FORCE
Mohammed,
please download the shapefile and kml from the links I have posted above. As mentioned before: it is not possible to cover your study area with S2 images from one single day.

What you posted are sentinel tile identifiers, not scene IDs. Again: check how orbits and S2 tiles overlay.

moh713...@gmail.com

unread,
Oct 27, 2020, 11:46:43 AM10/27/20
to FORCE
hi
so sorry  for this confusion.
some thing has gone wrong with me !!
I am reprocess the images again.
so sorry again for this confusion
Reply all
Reply to author
Forward
0 new messages