How to use the nansat.tools.gdal.GCP function in nansat

gcpLineName = '%s_%03d' % (gcpName, n)
                gcpString += geoMetadata[gcpLineName]
            # convert strings to floats
            gcpString = gcpString.strip().replace(' ', '')
            gcpValues = []
            # append all gcps from string
            for x in gcpString.split('|'):
                if len(x) > 0:
                    gcpValues.append(float(x))
            # gcpValues = [float(x) for x in gcpString.strip().split('|')]
            gcpAllValues.append(gcpValues)

        # create list of GDAL GCPs
        gcps = []
        for i in range(0, len(gcpAllValues[0])):
            gcps.append(gdal.GCP(gcpAllValues[2][i], gcpAllValues[3][i], 0,
                                 gcpAllValues[0][i], gcpAllValues[1][i]))

        return gcps

gcpLineName = '%s_%03d' % (gcpName, n)
                gcpString += geoMetadata[gcpLineName]
            # convert strings to floats
            gcpString = gcpString.strip().replace(' ', '')
            gcpValues = []
            # append all gcps from string
            for x in gcpString.split('|'):
                if len(x) > 0:
                    gcpValues.append(float(x))
            #gcpValues = [float(x) for x in gcpString.strip().split('|')]
            gcpAllValues.append(gcpValues)

        # create list of GDAL GCPs
        gcps = []
        for i in range(0, len(gcpAllValues[0])):
            gcps.append(gdal.GCP(gcpAllValues[2][i], gcpAllValues[3][i], 0,
                                 gcpAllValues[0][i], gcpAllValues[1][i]))

        if len(gcps) > 0:
            # get GCP projection and repare
            projection = self.repare_projection(geoMetadata.
                                                get('GCPProjection', ''))
            # add GCPs to dataset
            self.dataset.SetGCPs(gcps, projection)
            self._remove_geotransform()

        return len(gcps)

""" Create GCPs from geolocation data

        Parameters
        ----------
        x, y, z, p, l
        N-D arrays with value of X, Y, Z, Pixel and Line coordinates.
        X and Y are typically lon, lat, Z - height.

        Returns
        -------
        gcps : list with GDAL GCPs

        """
        gcps = []
        for xi, yi, zi, pi, li in zip(x.flat, y.flat, z.flat, p.flat, l.flat):
            gcps.append(gdal.GCP(xi, yi, zi, pi, li))
        return gcps

title,
                          latitude.shape[0], latitude.shape[1],
                          GCP_COUNT, step0, step1)

        # generate list of GCPs
        dx = .5
        dy = .5
        gcps = []
        k = 0
        for i0 in range(0, latitude.shape[0], step0):
            for i1 in range(0, latitude.shape[1], step1):
                # create GCP with X,Y,pixel,line from lat/lon matrices
                lon = float(longitude[i0, i1])
                lat = float(latitude[i0, i1])
                if (lon >= -180 and lon <= 180 and lat >= -90 and lat <= 90):
                    gcp = gdal.GCP(lon, lat, 0, i1 * pixelStep + dx,
                                   i0 * lineStep + dy)
                    self.logger.debug('%d %d %d %f %f',
                                      k, gcp.GCPPixel, gcp.GCPLine,
                                      gcp.GCPX, gcp.GCPY)
                    gcps.append(gcp)
                    k += 1

        # append GCPs and lat/lon projection to the vsiDataset
        self.dataset.SetGCPs(gcps, NSR().wkt)

GCP_COUNT, step0, step1)

        # generate list of GCPs
        dx = .5
        dy = .5
        gcps = []
        k = 0
        center_lon = 0
        center_lat = 0
        for i0 in range(0, latitude.shape[0], step0):
            for i1 in range(0, latitude.shape[1], step1):
                # create GCP with X,Y,pixel,line from lat/lon matrices
                lon = float(longitude[i0, i1])
                lat = float(latitude[i0, i1])
                if (lon >= -180 and lon <= 180 and lat >= -90 and lat <= 90):
                    gcp = gdal.GCP(lon, lat, 0, i1 * pixelStep + dx,
                                   i0 * lineStep + dy)
                    self.logger.debug('%d %d %d %f %f',
                                      k, gcp.GCPPixel, gcp.GCPLine,
                                      gcp.GCPX, gcp.GCPY)
                    gcps.append(gcp)
                    center_lon += gcp.GCPX
                    center_lat += gcp.GCPY
                    k += 1


        # append GCPs and lat/lon projection to the vsiDataset
        self.dataset.SetGCPs(gcps, NSR().wkt)
        self._remove_geolocation()

        # reproject GCPs
        center_lon /= k

float(top_left_lat),
                       0,
                       0,
                       subDataset.RasterYSize)
        gcps.append(gcp)
        #self.logger.debug('%d %d %d %f %f', 2, gcp.GCPPixel, gcp.GCPLine,
        #                  gcp.GCPX, gcp.GCPY)
        gcp = gdal.GCP(float(top_right_lon),
                       float(top_right_lat),
                       0,
                       subDataset.RasterXSize,
                       subDataset.RasterYSize)
        gcps.append(gcp)
        #self.logger.debug('%d %d %d %f %f', 3, gcp.GCPPixel, gcp.GCPLine,
        #                  gcp.GCPX, gcp.GCPY)
        gcp = gdal.GCP(float(center_lon),
                       float(center_lat),
                       0,
                       int(np.round(subDataset.RasterXSize/2.)),
                       int(round(subDataset.RasterYSize/2.)))
        gcps.append(gcp)
        #self.logger.debug('%d %d %d %f %f', 4, gcp.GCPPixel, gcp.GCPLine,
        #                  gcp.GCPX, gcp.GCPY)

        # append GCPs and lat/lon projection to the vsiDataset
        latlongSRS = osr.SpatialReference()
        latlongSRS.ImportFromProj4("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
        latlongSRSWKT = latlongSRS.ExportToWkt()

        # create empty VRT dataset with geolocation only
        VRT.__init__(self,
                     srcRasterXSize=subDataset.RasterXSize,

numOfGCPs = i

            if numOfGCPs < 100:
                # create new 100 GPCs (10 x 10 regular matrix)
                pixArray = []
                linArray = []
                for newPix in np.r_[0:xSize:10j]:
                    for newLin in np.r_[0:ySize:10j]:
                        pixArray.append(newPix + xOff)
                        linArray.append(newLin + yOff)

                lonArray, latArray = self.vrt.transform_points(pixArray,
                                                               linArray)

                for i in range(len(lonArray)):
                    dstGCPs.append(gdal.GCP(lonArray[i], latArray[i], 0,
                                            pixArray[i] - xOff,
                                            linArray[i] - yOff,
                                            '', str(numOfGCPs+i+1)))

            # set new GCPss
            self.vrt.dataset.SetGCPs(dstGCPs, NSR().wkt)
            # reproject new GCPs to the SRS of original GCPs
            nsr = NSR(self.vrt.dataset.GetGCPProjection())
            self.vrt.reproject_GCPs(nsr)
            # remove geotranform which was automatically added
            self.vrt._remove_geotransform()
        else:
            # shift upper left corner coordinates
            geoTransfrom = self.vrt.dataset.GetGeoTransform()
            geoTransfrom = map(float, geoTransfrom)
            geoTransfrom[0] += geoTransfrom[1] * xOff

lonSubdataset = [subdatasetName[0]
                         for subdatasetName in gdalDataset.GetSubDatasets()
                         if 'Longitude' in subdatasetName[1]][0]
        latSubdataset = [subdatasetName[0]
                         for subdatasetName in gdalDataset.GetSubDatasets()
                         if 'Latitude' in subdatasetName[1]][0]
        lons = gdal.Open(lonSubdataset).ReadAsArray()
        lats = gdal.Open(latSubdataset).ReadAsArray()
        gcps = []
        rows = range(0, lons.shape[0], lons.shape[0]/GCP_COUNT)
        cols = range(0, lons.shape[1], lons.shape[1]/GCP_COUNT)
        factor = self.dataset.RasterYSize / lons.shape[0]
        for r in rows:
            for c in cols:
                gcps.append(gdal.GCP(float(lons[r,c]),
                                     float(lats[r,c]), 0,
                                     factor*c+0.5,
                                     factor*r+0.5))
        self.dataset.SetGCPs(gcps, self.dataset.GetGCPProjection())
        self.tps = True