GeoscienceAustralia / uncover-ml

Machine Learning system for Geoscience Australia uncover project
Apache License 2.0
30 stars 20 forks source link

Raster prediction output Affine transform y is positive, not negative #125

Open bluetyson opened 3 years ago

bluetyson commented 3 years ago

here's an output from the other day:-

with rasterio.open(r'F:\'RFNoNaN074_prediction.tif') as src:
    print(src.meta)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -1.0000000200408773e+20, 'width': 1322, 'height': 1596, 'count': 1, 'crs': CRS.from_epsg(3107), 'transform': Affine(1000.0, 0.0, 339300.0,
       0.0, 1000.0, 1136900.0)}

with rasterio.open(r'F:\SA_GRAV.tif') as src2:
    print(src2.meta)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -99999.0, 'width': 1322, 'height': 1596, 'count': 1, 'crs': CRS.from_epsg(3107), 'transform': Affine(1000.0, 0.0, 339300.0,
       0.0, -1000.0, 2732900.0)}

Noticed this when an output was read into ER Mapper. Thought maybe an oddball project thing as fine in QGIS and Arc.

However, Raster in R - same problem, so definitely the raster.

Then I checked, as above.

bluetyson commented 3 years ago

ymax and ymin around the wrong way in geoio.py?

self.A, _, _ = image.bbox2affine(bbox[1, 0], bbox[0, 0],
                                         bbox[0, 1], bbox[1, 1],
                                         shape[0], shape[1])

change that though will leave flipped data I think?

bluetyson commented 3 years ago

Might have to do something like this then?

self.sub_starts = np.flip(self.sub_starts)

and

for i, f in enumerate(self.files):
                        #print(i, f)
                        f.write(np.rot90(data[i:i+1],2), window=window)
RichardScottOZ commented 3 years ago

def __init__(self, shape, bbox, crs, n_subchunks, outpath,
                 outbands, band_tags=None, independent=False, **kwargs):`
        """
        pass in additional geotif write options in kwargs`
        """
       # affine`
       #print("bbox: ", bbox)`
       #print("bbox00: ", bbox[0,0])`
        #print("bbox01: ", bbox[0,1])`
        #print("bbox10: ", bbox[1,0])`
        #print("bbox11: ", bbox[1,1])`
        #xmax, xmin, ymax, ymin`
        #self.A, _, _ = image.bbox2affine(bbox[1, 0], bbox[0, 0],`
        #                                 bbox[0, 1], bbox[1, 1],`
        #                                 shape[0], shape[1])`
        self.A, _, _ = image.bbox2affine(bbox[1, 0], bbox[0, 0],  bbox[1, 1], bbox[0, 1],   shape[0], shape[1])`
        #print("imagebbox2affine: ",image.bbox2affine(bbox[1, 0], bbox[0, 0],`
        #                                 bbox[0, 1], bbox[1, 1],`
        #                                 shape[0], shape[1]))`

        #print("ibbabb: ",bbox[1, 0], bbox[0, 0], bbox[0, 1], bbox[1, 1], shape[0], shape[1])`

        self.shape = shape
        self.outbands = outbands
        self.bbox = bbox
        self.outpath = outpath
        self.n_subchunks = n_subchunks
        self.independent = independent  # mpi control
        self.sub_starts = [k[0] for k in np.array_split(
                           np.arange(self.shape[1]),
                           mpiops.chunks * self.n_subchunks)]
        #print("sub_starts",self.sub_starts)
        self.sub_starts = np.flip(self.sub_starts)  #flipping for new transform fix
        #print("new sub_starts",self.sub_starts)

        # file tags don't have spaces
        if band_tags:
            if self.outbands > len(band_tags):
                _logger.warning(f"Specified more outbands ({self.outbands}) than there are "
                                f"prediction tags available ({len(band_tags)}). "
                                f"Limiting outbands to number of prediction tags.")
                self.outbands = len(band_tags)
            file_tags = ["_".join(k.lower().split()) for k in band_tags]
        else:
            file_tags = [str(k) for k in range(self.outbands)]
            band_tags = file_tags

        files = []
        file_names = []

        if mpiops.chunk_index == 0:
            for band in range(self.outbands):
                output_filename = self.outpath.format(file_tags[band])
                f = rasterio.open(output_filename, 'w', driver='GTiff',
                                  width=self.shape[0], height=self.shape[1],
                                  dtype=np.float32, count=1,
                                  crs=crs,
                                  transform=self.A,
                                  nodata=self.nodata_value,
                                  **kwargs
                                  )
                f.update_tags(1, image_type=band_tags[band])
                files.append(f)
                file_names.append(output_filename)

        if independent:
            self.files = files
        else:
            if mpiops.chunk_index == 0:
                # create a file for each band
                self.files = files
                self.file_names = file_names
            else:
                self.file_names = []

            self.file_names = mpiops.comm.bcast(self.file_names, root=0)
RichardScottOZ commented 2 years ago

Transform is still the same incorrect on latest version.