Hello everyone !
We are french students in remote sensing and for a project we try to classify ice, open water and ground in the Arctic region (North of Svalbard Islands).
We use Sentinel 1 GRD EW data. We are trying to automate pre-processing. Before we would like to make a subset using only one coordinate point. With this point, we add an area of dozens of kilometers (the point is in the middle). So, we manually enter geocoordinates in the subset tool.
When we launch the Python script, the result is unusable and we think it’s maybe because of the subset. In fact, we tried to do the same process without the subset and it works perfectly. We also tried to draw manually the subset with the blue rectangle and it worked.
So, we were wondering, with which EPSG does the subset tool work? Is it WGS 84?
Here are our coordinates:
longitude_min: 10.912
longitude_max: 7.610
latitude_min: 77.904
latitude_max: 78.625
And here is our script:
import os
import glob
import subprocess
from sentinelsat import SentinelAPI
from pyproj import Proj, transform
from shapely.geometry import Point
with open("D:\\Projet_pro\\coord.txt") as coord_vincent:
# Extraction of the point's coordinates
coord = coord_vincent.read()
lat, lon = coord.split(",")
print(f"Longitude : {lon}° - Latitude : {lat}°")
# Conversion in EPSG Polar stereographic
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:3995')
x, y = transform(inProj, outProj, lon, lat)
print(f"x : {x} m - y : {y} m")
#Add a perimeter in meters
xmin = x-30000
ymin = y-15000
xmax = x+30000
ymax = y+75000
#Reprojection in EPSG WGS84 to use in SNAP
inProj = Proj(init='epsg:3995')
outProj = Proj(init='epsg:4326')
lonmax, latmin = transform(inProj, outProj, xmin, ymin)
lonmin, latmax = transform(inProj, outProj, xmax, ymax)
print(lonmin,lonmax,latmin,latmax)
rep="D:\\Projet_pro" #Project directoryretrouver les graphs.xml, les répertoires et le geojson
modelgraph="ModelProcessSNAP.xml" #Name of the initial Graph Builder
graph="myGraph.xml" #Name of the complete graph written for each image
snap="C:\\Program Files\\snap\\bin\\gpt" #SNAP directory
repimg = "D:\\Projet_pro\\Download" #File for downloaded images
repoutput = "D:\\Projet_pro\\Output" #File for preprocessed images
#Traitements
files = glob.glob(os.path.join(repimg,'*.zip')) # Listing of each .zip in this repertory
for f in files:
# Transform the modelgraph.xml in a complete Graph Builder
foutname = os.path.join(repoutput, os.path.basename(f)[:-4]+"_output.dim")
infile = open(os.path.join(rep,modelgraph), 'r')
tempFile = open(os.path.join(rep,graph), 'w')
for line in infile:
if "FILEIN" in line :
line = line.replace( "FILEIN", f )#<- HERE, we replace "FILEIN" by the name of the file to treat
elif "FILEOUT" in line :
line = line.replace( "FILEOUT", foutname)
elif "COORDIN" in line:
scoord = f"{lonmin} {latmax}, {lonmax} {latmax}, {lonmax} {latmin}, {lonmin} {latmin}, {lonmin} {latmax}, {lonmin} {latmax} "
line = line.replace("COORDIN", scoord)
tempFile.write(line)
tempFile.close()
infile.close()
#Execute the graph in SNAP
print("### Execution SNAP pour fichier '" + f + "' en cours ...")
p = subprocess.Popen([snap, os.path.join(rep,graph)], stdout=subprocess.PIPE)
print("------ Trace SNAP -------")
out = p.stdout.read()
print(out)
print("---------------------------")
Thanks for you help.