About MSI's output results

Post Reply
Posts: 1
Joined: Sun Oct 23, 2022 5:51 pm
company / institution: Tsukuba University
Location: Japan

About MSI's output results

Post by MAILISU »

Dear Steinmetz,

Thank you for Polymer.
I'm working on a matching of two MSI images(10m) after polymer processing.
I stitching the two images in the code with the same rows and columns(3300*3600).
However, When I stitch the two images together there is a large shift in the area where they should overlap.
I checked the Geotransform and projection(both files have the same projections):
Geotransform 1: (140.1548560229316, 0.00011158006136450954, 0.0, 36.14402615549123, 0.0, -9.077539856260909e-05)
Geotransform 2: (140.15259911600464, 0.00011188156175087727, 0.0, 36.35336131677865, 0.0, -9.077718296325087e-05)
Latitude and longitude of the file:
file1:140.1548560229316 35.84446734023462 140.55654424384383 36.14402615549123
file2:140.15259911600464 36.053796612999925 140.5553727383078 36.35336131677865
Images of the same size differ in longitude by 0.001

and I checked the latitude and longitude of the overlapping features and the difference in longitude of the pixels at the same location is almost 0.002 degrees and the difference in latitude is 0.001 degrees.
I also processed other images from the same location with different dates, all producing the same displacement.
Is this due to a problem with polymer's original code that causes the latitude and longitude of the two images to be inconsistent?

Thank you very much!

(I've attached the code I used to convert the nc file to tiff format, please correct me if this code is causing problems with the images not overlapping.)

def Convert_netCDF_to_GTiff_polymer(netCDF_File_List, overwrite=False):
GTiff_File_List = []
para = []
for i, ifile in enumerate(netCDF_File_List):
tic = datetime.now()
print(f'Processing file no.{i + 1:02d} (/{len(Sat_Data_info)})...')
### reading info from netCDF file
f = Dataset(ifile, 'r')
lon, lat = f.variables['longitude'][:], f.variables['latitude'][:]
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
print(xmin, ymin, xmax, ymax)
epsg = get_UTM_EPSGcode(lon, lat)
row, col = lon.shape
xRes, yRes = (xmax - xmin) / col, (ymax - ymin) / row
geotransform = (xmin, xRes, 0, ymax, 0, -yRes)#
para = [ip for ip in f.variables if ('lat' not in ip) & ('lon' not in ip)]
# ptype = [f.variables[ip].dtype for ip in para]
fillvalue = [f.variables[ip]._FillValue for ip in para]
FillValue = -32767.
nband = len(para)

### writting to GeoTiff file
ofile = ifile.replace('.nc', '.tiff')
if (os.path.exists(ofile) == True) & (overwrite == False):
print('GTiff file exists. SKIP\n')
dst = gdal.GetDriverByName('GTiff').Create(ofile, col, row, nband, gdal.GDT_Float32)
dst.SetGeoTransform(geotransform) # specify coords
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg) # WGS84 lat/long
dst.SetProjection(srs.ExportToWkt()) # export coords to file
for ii, (ip, ifv) in enumerate(zip(para, fillvalue)):
print(f'writing {ip} to Band{ii + 1:02d}...')
# The TIFFTAG_GDAL_NODATA only support one value per dataset, so -32767 is assigned to all bands
data = np.array(f.variables[ip][:]).astype(np.float32)
data = np.where(data == ifv, FillValue, data)
dst.GetRasterBand(ii + 1).WriteArray(data)
dst.GetRasterBand(ii + 1).SetNoDataValue(FillValue)
dst = None
toc = datetime.now()
print(f'Elapsed time is {(toc - tic).total_seconds():0.2f} seconds\n')

return para
Post Reply