Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

AttributeError: 'NoneType' object has no attribute 'to_epsg' #181

Open
weiji14 opened this issue Sep 23, 2022 · 9 comments · May be fixed by #182
Open

AttributeError: 'NoneType' object has no attribute 'to_epsg' #181

weiji14 opened this issue Sep 23, 2022 · 9 comments · May be fixed by #182

Comments

@weiji14
Copy link

weiji14 commented Sep 23, 2022 •
edited
Loading

Hi there, just encountering some issues while trying to stack Sentinel-1 GRD data from Planetary Computer at weiji14/zen3geo#62, and I'm wondering if it's an issue on the STAC metadata side (c.f. #152), or if it's something to fix in stackstac.

Here's a MCWE, using pystac=1.4.0, planetary-computer=0.4.7, stackstac=0.4.3:

import stackstac
import pystac
import planetary_computer
import xarray as xr

# %%
item_urls: list = [
    "https://planetarycomputer.microsoft.com/api/stac/v1/collections/sentinel-1-grd/items/S1A_IW_GRDH_1SDV_20220320T230514_20220320T230548_042411_050E99",
    "https://planetarycomputer.microsoft.com/api/stac/v1/collections/sentinel-1-grd/items/S1A_IW_GRDH_1SDV_20220308T230513_20220308T230548_042236_0508AF",
]

# Load the individual item metadata and sign the assets
items: list[pystac.Item] = [pystac.Item.from_file(item_url) for item_url in item_urls]
signed_items: list[pystac.Item] = [planetary_computer.sign(item) for item in items]

# Stack Sentinel-1 GRD files
dataarray: xr.DataArray = stackstac.stack(
    items=signed_items,
    epsg=32647,
    resolution=30,
    bounds_latlon=[99.933681, -0.009951, 100.065765, 0.147054],  # W, S, E, N
)
assert dataarray.crs == "epsg:32647"
print(dataarray)

The output xarray repr looks ok like so. Dimensions are (time:2, band:2, y:579, x:491).

<xarray.DataArray 'stackstac-950602eb423dd1d439106f6794699f05' (time: 2,
                                                                band: 2,
                                                                y: 579, x: 491)>
dask.array<fetch_raster_window, shape=(2, 2, 579, 491), dtype=float64, chunksize=(1, 1, 579, 491), chunktype=numpy.ndarray>
Coordinates: (12/37)
  * time                                   (time) datetime64[ns] 2022-03-08T2...
    id                                     (time) <U62 'S1A_IW_GRDH_1SDV_2022...
  * band                                   (band) <U2 'vh' 'vv'
  * x                                      (x) float64 6.039e+05 ... 6.186e+05
  * y                                      (y) float64 1.626e+04 ... -1.08e+03
    end_datetime                           (time) <U32 '2022-03-08 23:05:48.2...
    ...                                     ...
    s1:resolution                          <U4 'high'
    s1:product_timeliness                  <U8 'Fast-24h'
    sat:relative_orbit                     int64 164
    description                            (band) <U145 'Amplitude of signal ...
    title                                  (band) <U41 'VH: vertical transmit...
    epsg                                   int64 32647
Attributes:
    spec:        RasterSpec(epsg=32647, bounds=(603870, -1110, 618600, 16260)...
    crs:         epsg:32647
    transform:   | 30.00, 0.00, 603870.00|\n| 0.00,-30.00, 16260.00|\n| 0.00,...
    resolution:  30

but when I try to run dataarray.compute(), this AttributeError: 'NoneType' object has no attribute 'to_epsg' message popped up. Here's the full traceback:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 dataarray.compute()

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/xarray/core/dataarray.py:993, in DataArray.compute(self, **kwargs)
    974 """Manually trigger loading of this array's data from disk or a
    975 remote source into memory and return a new array. The original is
    976 left unaltered.
   (...)
    990 dask.compute
    991 """
    992 new = self.copy(deep=False)
--> 993 return new.load(**kwargs)

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/xarray/core/dataarray.py:967, in DataArray.load(self, **kwargs)
    949 def load(self: T_DataArray, **kwargs) -> T_DataArray:
    950     """Manually trigger loading of this array's data from disk or a
    951     remote source into memory and return this array.
    952 
   (...)
    965     dask.compute
    966     """
--> 967     ds = self._to_temp_dataset().load(**kwargs)
    968     new = self._from_temp_dataset(ds)
    969     self._variable = new._variable

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/xarray/core/dataset.py:733, in Dataset.load(self, **kwargs)
    730 import dask.array as da
    732 # evaluate all the dask arrays simultaneously
--> 733 evaluated_data = da.compute(*lazy_data.values(), **kwargs)
    735 for k, data in zip(lazy_data, evaluated_data):
    736     self.variables[k].data = data

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/base.py:602, in compute(traverse, optimize_graph, scheduler, get, *args, **kwargs)
    599     keys.append(x.__dask_keys__())
    600     postcomputes.append(x.__dask_postcompute__())
--> 602 results = schedule(dsk, keys, **kwargs)
    603 return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/threaded.py:89, in get(dsk, keys, cache, num_workers, pool, **kwargs)
     86     elif isinstance(pool, multiprocessing.pool.Pool):
     87         pool = MultiprocessingPoolExecutor(pool)
---> 89 results = get_async(
     90     pool.submit,
     91     pool._max_workers,
     92     dsk,
     93     keys,
     94     cache=cache,
     95     get_id=_thread_get_id,
     96     pack_exception=pack_exception,
     97     **kwargs,
     98 )
    100 # Cleanup pools associated to dead threads
    101 with pools_lock:

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/local.py:511, in get_async(submit, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, chunksize, **kwargs)
    509         _execute_task(task, data)  # Re-execute locally
    510     else:
--> 511         raise_exception(exc, tb)
    512 res, worker_id = loads(res_info)
    513 state["cache"][key] = res

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/local.py:319, in reraise(exc, tb)
    317 if exc.__traceback__ is not tb:
    318     raise exc.with_traceback(tb)
--> 319 raise exc

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/local.py:224, in execute_task(key, task_info, dumps, loads, get_id, pack_exception)
    222 try:
    223     task, data = loads(task_info)
--> 224     result = _execute_task(task, data)
    225     id = get_id()
    226     result = dumps((result, id))

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/core.py:119, in _execute_task(arg, cache, dsk)
    115     func, args = arg[0], arg[1:]
    116     # Note: Don't assign the subtask results to a variable. numpy detects
    117     # temporaries by their reference count and can execute certain
    118     # operations in-place.
--> 119     return func(*(_execute_task(a, cache) for a in args))
    120 elif not ishashable(arg):
    121     return arg

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/optimization.py:990, in SubgraphCallable.__call__(self, *args)
    988 if not len(args) == len(self.inkeys):
    989     raise ValueError("Expected %d args, got %d" % (len(self.inkeys), len(args)))
--> 990 return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/core.py:149, in get(dsk, out, cache)
    147 for key in toposort(dsk):
    148     task = dsk[key]
--> 149     result = _execute_task(task, cache)
    150     cache[key] = result
    151 result = _execute_task(out, cache)

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/dask/core.py:119, in _execute_task(arg, cache, dsk)
    115     func, args = arg[0], arg[1:]
    116     # Note: Don't assign the subtask results to a variable. numpy detects
    117     # temporaries by their reference count and can execute certain
    118     # operations in-place.
--> 119     return func(*(_execute_task(a, cache) for a in args))
    120 elif not ishashable(arg):
    121     return arg

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/stackstac/to_dask.py:185, in fetch_raster_window(reader_table, slices, dtype, fill_value)
    178 # Only read if the window we're fetching actually overlaps with the asset
    179 if windows.intersect(current_window, asset_window):
    180     # NOTE: when there are multiple assets, we _could_ parallelize these reads with our own threadpool.
    181     # However, that would probably increase memory usage, since the internal, thread-local GDAL datasets
    182     # would end up copied to even more threads.
    183 
    184     # TODO when the Reader won't be rescaling, support passing `output` to avoid the copy?
--> 185     data = reader.read(current_window)
    187     if all_empty:
    188         # Turn `output` from a broadcast-trick array to a real array, so it's writeable
    189         if (
    190             np.isnan(data)
    191             if np.isnan(fill_value)
    192             else np.equal(data, fill_value)
    193         ).all():
    194             # Unless the data we just read is all empty anyway

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/stackstac/rio_reader.py:385, in AutoParallelRioReader.read(self, window, **kwargs)
    384 def read(self, window: Window, **kwargs) -> np.ndarray:
--> 385     reader = self.dataset
    386     try:
    387         result = reader.read(
    388             window=window,
    389             masked=True,
   (...)
    392             **kwargs,
    393         )

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/stackstac/rio_reader.py:381, in AutoParallelRioReader.dataset(self)
    379 with self._dataset_lock:
    380     if self._dataset is None:
--> 381         self._dataset = self._open()
    382     return self._dataset

File ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/stackstac/rio_reader.py:348, in AutoParallelRioReader._open(self)
    340     raise RuntimeError(
    341         f"Assets must have exactly 1 band, but file {self.url!r} has {ds.count}. "
    342         "We can't currently handle multi-band rasters (each band has to be "
    343         "a separate STAC asset), so you'll need to exclude this asset from your analysis."
    344     )
    346 # Only make a VRT if the dataset doesn't match the spatial spec we want
    347 if self.spec.vrt_params != {
--> 348     "crs": ds.crs.to_epsg(),
    349     "transform": ds.transform,
    350     "height": ds.height,
    351     "width": ds.width,
    352 }:
    353     with self.gdal_env.open_vrt:
    354         vrt = WarpedVRT(
    355             ds,
    356             sharing=False,
    357             resampling=self.resampling,
    358             **self.spec.vrt_params,
    359         )

AttributeError: 'NoneType' object has no attribute 'to_epsg'

So I verified that the stacked Sentinel-1 dataarray does have a crs attribute like 'epsg:32647', but it seems like ds is something else? I did try to step through the code at

stackstac/stackstac/rio_reader.py

Lines 347 to 362 in 9106708

if self.spec.vrt_params != {
"crs": ds.crs.to_epsg(),
"transform": ds.transform,
"height": ds.height,
"width": ds.width,
}:
with self.gdal_env.open_vrt:
vrt = WarpedVRT(
ds,
sharing=False,
resampling=self.resampling,
**self.spec.vrt_params,
)
else:
logger.info(f"Skipping VRT for {self.url!r}")
vrt = None

but am having trouble understaing what ds is representing, or how ds.crs can be set to a proper coordinate reference system value. My guess that the crs needs to be set in the STAC Item metadata? Or perhaps the if-statement needs to be revised.

Oh, and just to make sure the Sentinel-1 GRD STAC Items are readable, I did try reading it using rioxarray=0.11.1:

import rioxarray

url: str = signed_items[1].assets["vv"].href
dataarray = rioxarray.open_rasterio(filename=url, overview_level=5)
dataarray = dataarray.compute()
print(dataarray)
# <xarray.DataArray (band: 1, y: 361, x: 396)>
# array([[[  0, 224, 259, ...,  45,   0,   0],
#         [  0, 243, 286, ...,  44,   0,   0],
#         [  0, 274, 248, ...,  43,   0,   0],
#         ...,
#         [  0,   0,   0, ...,  34,  36,  36],
#         [  0,   0,   0, ...,  33,  35,  34],
#         [  0,   0,   0, ...,  17,  17,  17]]], dtype=uint16)
# Coordinates:
#   * band         (band) int64 1
#     spatial_ref  int64 0
# Dimensions without coordinates: y, x
# Attributes:
#     _FillValue:    0.0
#     scale_factor:  1.0
#     add_offset:    0.0

Here's the output from rioxarray.show_version() for completeness to show my GDAL and other geo-library versions:

rioxarray (0.11.1) deps:
  rasterio: 1.3.0
    xarray: 2022.3.0
      GDAL: 3.5.0
      GEOS: 3.10.2
      PROJ: 9.0.0
 PROJ DATA: ~/mambaforge/envs/zen3geo/lib/python3.10/site-packages/rasterio/proj_data
 GDAL DATA: None

Other python deps:
     scipy: 1.9.0
    pyproj: 3.3.1

System:
    python: 3.10.6 | packaged by conda-forge | (main, Aug 22 2022, 20:36:39) [GCC 10.4.0]
executable: ~/mambaforge/envs/zen3geo/bin/python
   machine: Linux-5.17.0-1016-oem-x86_64-with-glibc2.35

P.S. Thanks for this amazing library, really like the design and it's been a pleasure to use 😄
@TomAugspurger
Copy link
Contributor

TomAugspurger commented Sep 23, 2022 •
edited
Loading

The root issue is that the rasterio doesn't read a CRS from the COG:

import pystac
import rasterio
import planetary_computer

item = planetary_computer.sign(pystac.read_file("https://planetarycomputer.microsoft.com/api/stac/v1/collections/sentinel-1-grd/items/S1A_IW_GRDH_1SDV_20220320T230514_20220320T230548_042411_050E99"))
rasterio.open(item.assets["vh"].href).crs  # None

I'll look into it. I think it was a deliberate choice because of something about the data.

FYI @weiji14, if the RTC product is OK for your application, it does have the CRS (and I think something about the RTC process is what makes having a CRS for it correct, but not for the GRD).

@lossyrob
Copy link

Sentinel-1 GRD does not have a referenceable CRS because it uses Ground Control Points (GCPs), which ties a subset of pixels to points on the surface of the Earth, but do not necessarily conform to any well-known projection. The GRD product is a precursor to the radiometric correction and terrain flattening that is required in order to map each pixel value to distinct points on the ground. The RTC product referenced above is the result of that processing, and it's recommended to use the RTC product when possible for applications that need georectified raster data.

@weiji14
Copy link
Author

weiji14 commented Sep 23, 2022

Oh yes, you're both right. The Sentinel-1 GRD x and y values (read from rioxarray) are actually just image coordinates and not geographical coordinates. Somehow I thought the GRD product was already georeferenced from the SLC (or semi-georeferenced), since it shows up nicely on Planetary Computer Explorer, but apparently not! This makes me wonder though, is there some affine transform in the STAC metadata I can use to convert the GRD image coordinates to some rough geographical coordinates?

FYI @weiji14, if the RTC product is OK for your application, it does have the CRS (and I think something about the RTC process is what makes having a CRS for it correct, but not for the GRD).

Yes, I'm aware of the RTC dataset but 1) it requires authentication and I haven't been able to get it outside of the Planetary Computer environment (even though I tried setting the PC_SDK_SUBSCRIPTION_KEY following https://planetarycomputer.microsoft.com/dataset/sentinel-1-rtc#Example-Notebook) and 2) I'm trying to get this to work on a notebook hosted publicly on readthedocs, and would rather not mess with secret auth tokens 🙂

@scottyhq
Copy link
Contributor

^ Agreed that RTC is generally the better choice for analysis! GRD may be okay for certain analyses though, especially in flat terrain.

@weiji14 GCPs to well-known CRS must make interpolation choices like what order polynomial, thin-plate-spline, etc. Which you can definitely do with rioxarray and rasterio>1.2 (see more detail in corteva/rioxarray#339 (comment) or basic behavior below):

# GRD in geographic coords using basic GCP polynomial:
with rasterio.open(grd_url) as src:
    print(src.crs) #None
    with rasterio.vrt.WarpedVRT(src) as vrt:
        print(vrt.crs) #EPSG:4326
        da_grd = rioxarray.open_rasterio(vrt)
        print(da_grd.rio.crs) #EPSG:4326

The current stackstack behavior does seem confusing. 1. stackstac.stack could instead fail with "GCPs or lack of source CRS not supported". Or 2. since WarpedVRT is used by stackstack perhaps the basic GCP transform illustrated above could be used by default...

@weiji14
Copy link
Author

weiji14 commented Sep 23, 2022 •
edited
Loading

@weiji14 GCPs to well-known CRS must make interpolation choices like what order polynomial, thin-plate-spline, etc. Which you can definitely do with rioxarray and rasterio>1.2 (see more detail in corteva/rioxarray#339 (comment) or basic behavior below):

Ah that's good to know (thanks for chiming in @scottyhq!). I see that the Sentinel-1 GRD does have GCPs mapping pixel row/col coordinates to geographic lon/lat coordinates.

with rasterio.open(grd_url) as src:
    print(src.gcps)
# ([GroundControlPoint(row=0.0, col=0.0, x=100.6583362878721, y=-0.010462088462186858, z=610.0002990607172, id='1', info=''), GroundControlPoint(row=0.0, col=1266.0, x=100.54377974847928, y=0.014154768329878284, z=840.0004153857008, id='2', info=''), GroundControlPoint(row=0.0, col=2532.0, x=100.42944973594832, y=0.03872256105283347, z=1070.0005301199853, id='3', info=''), GroundControlPoint(row=0.0, col=3798.0, x=100.31638546987982, y=0.06301688895194679, z=1220.0006033526734, id='4', info=''), GroundControlPoint(row=0.0, col=5064.0, x=100.2083214744807, y=0.08623172797820164, z=990.000489811413, id='5', info=''), GroundControlPoint(row=0.0, col=6330.0, x=100.09913352366883, y=0.10968833549442128, z=835.0004128646106, id='6', info=''), GroundControlPoint(row=0.0, col=7596.0, x=99.99346436439382, y=0.13238464099033495, z=380.0001940652728, id='7', info=''), GroundControlPoint(row=0.0, col=8862.0, x=99.88320268867845, y=0.15607110859024764, z=288.0001501413062, id='8', info=''), GroundControlPoint(row=0.0, col=10128.0, x=99.77169560224651, y=0.18002542944469227, z=300.0001548547298, id='9', info=''), GroundControlPoint(row=0.0, col=11394.0, x=99.66267664792453, y=0.2034415498912254, z=90.00005884934217, id='10', info=''), GroundControlPoint(row=0.0, col=12660.0, x=99.55126217744925, y=0.2273736206317498, z=90.00005904398859, id='11', info=''), GroundControlPoint(row=0.0, col=13926.0, x=99.43954055413064, y=0.2513707084396476, z=120.00007176958025, id='12', info=''), GroundControlPoint(row=0.0, col=15192.0, x=99.32876925137741, y=0.2751612325253469, z=60.00004603154957, id='13', info=''), GroundControlPoint(row=0.0, col=16458.0, x=99.2175234131717, y=0.29905283817730455, z=45.000039683654904, id='14', info=''), GroundControlPoint(row=0.0, col=17724.0, x=99.10657098908595, y=0.3228794966940388, z=2.1921470761299133e-05, id='15', info=''), GroundControlPoint(row=0.0, col=18990.0, x=98.99517536377506, y=0.3468001670665459, z=2.2017396986484528e-05, id='16', info=''), GroundControlPoint(row=0.0, col=20256.0, x=98.88378693029364, y=0.3707175628989745, z=2.19075009226799e-05, id='17', info=''), GroundControlPoint(row=0.0, col=21522.0, x=98.77240514147805, y=0.39463170333105724, z=2.1588988602161407e-05, id='18', info=''), GroundControlPoint(row=0.0, col=22788.0, x=98.66102949350878, y=0.4185425982084825, z=2.106931060552597e-05, id='19', info=''), GroundControlPoint(row=0.0, col=24054.0, x=98.54965950835697, y=0.4424501874633949, z=2.0342878997325897e-05, id='20', info=''), GroundControlPoint(row=0.0, col=25312.0, x=98.43899847716469, y=0.4662034713489156, z=1.942366361618042e-05, id='21', info=''), GroundControlPoint(row=2015.0, col=0.0, x=100.61914770478661, y=-0.1917016517096885, z=685.0004720995203, id='22', info=''), GroundControlPoint(row=2015.0, col=1266.0, x=100.5052152951129, y=-0.16719713289676746, z=885.000634457916, id='23', info=''), GroundControlPoint(row=2015.0, col=2532.0, x=100.39606630627804, y=-0.1437250901399487, z=750.0005434993654, id='24', info=''), GroundControlPoint(row=2015.0, col=3798.0, x=100.28845063476182, y=-0.12058403464251062, z=490.0003536650911, id='25', info=''), GroundControlPoint(row=2015.0, col=5064.0, x=100.18003331396312, y=-0.09726981993924248, z=275.0001949155703, id='26', info=''), GroundControlPoint(row=2015.0, col=6330.0, x=100.06812965197365, y=-0.0732025683287825, z=325.00023950915784, id='27', info=''), GroundControlPoint(row=2015.0, col=7596.0, x=99.95881760707738, y=-0.04969527375643036, z=165.00012012012303, id='28', info=''), GroundControlPoint(row=2015.0, col=8862.0, x=99.84780316064925, y=-0.025820310072807804, z=140.00010506808758, id='29', info=''), GroundControlPoint(row=2015.0, col=10128.0, x=99.73809599245253, y=-0.002228060531357331, z=-3.594905138015747e-07, id='30', info=''), GroundControlPoint(row=2015.0, col=11394.0, x=99.62674414093905, y=0.02171923193697134, z=3.3928081393241882e-06, id='31', info=''), GroundControlPoint(row=2015.0, col=12660.0, x=99.51539927131127, y=0.04566457095521927, z=6.925314664840698e-06, id='32', info=''), GroundControlPoint(row=2015.0, col=13926.0, x=99.40406083979661, y=0.06960810376542143, z=1.023896038532257e-05, id='33', info=''), GroundControlPoint(row=2015.0, col=15192.0, x=99.29272829048953, y=0.09354966276229774, z=1.3340264558792114e-05, id='34', info=''), GroundControlPoint(row=2015.0, col=16458.0, x=99.18140118807645, y=0.11748942644946991, z=1.622457057237625e-05, id='35', info=''), GroundControlPoint(row=2015.0, col=17724.0, x=99.07007907660423, y=0.14142727524580836, z=1.8890947103500366e-05, id='36', info=''), GroundControlPoint(row=2015.0, col=18990.0, x=98.95876155885787, y=0.16536319104837266, z=2.1342188119888306e-05, id='37', info=''), GroundControlPoint(row=2015.0, col=20256.0, x=98.84744828341816, y=0.18929721531514945, z=2.3577362298965454e-05, id='38', info=''), GroundControlPoint(row=2015.0, col=21522.0, x=98.7361389142482, y=0.21322932175103487, z=2.5598332285881042e-05, id='39', info=''), GroundControlPoint(row=2015.0, col=22788.0, x=98.62483314067364, y=0.23715947862279507, z=2.7406960725784302e-05, id='40', info=''), GroundControlPoint(row=2015.0, col=24054.0, x=98.51353066169835, y=0.26108758773617174, z=2.9003247618675232e-05, id='41', info=''), GroundControlPoint(row=2015.0, col=25312.0, x=98.40293452505902, y=0.2848624731448204, z=3.0377879738807678e-05, id='42', info=''), GroundControlPoint(row=4030.0, col=0.0, x=100.58213958969631, y=-0.3734153426767325, z=610.0005085738376, id='43', info=''), GroundControlPoint(row=4030.0, col=1266.0, x=100.46499377050984, y=-0.3481943796104414, z=1035.0009589577094, id='44', info=''), GroundControlPoint(row=4030.0, col=2532.0, x=100.35339555303858, y=-0.3241710127649846, z=1085.001035194844, id='45', info=''), GroundControlPoint(row=4030.0, col=3798.0, x=100.24348366517556, y=-0.30051122323295454, z=1010.0009805262089, id='46', info=''), GroundControlPoint(row=4030.0, col=5064.0, x=100.14167858012317, y=-0.27860361748731244, z=300.00024903565645, id='47', info=''), GroundControlPoint(row=4030.0, col=6330.0, x=100.03211925165611, y=-0.25501881572586793, z=165.0001152595505, id='48', info=''), GroundControlPoint(row=4030.0, col=7596.0, x=99.9224628764349, y=-0.23141238608885661, z=29.99997778236866, id='49', info=''), GroundControlPoint(row=4030.0, col=8862.0, x=99.81147709937484, y=-0.20751787007898984, z=-4.809629172086716e-05, id='50', info=''), GroundControlPoint(row=4030.0, col=10128.0, x=99.70013438810074, y=-0.18354564690562788, z=-4.1690655052661896e-05, id='51', info=''), GroundControlPoint(row=4030.0, col=11394.0, x=99.58879894560893, y=-0.15957452849885406, z=-3.551039844751358e-05, id='52', info=''), GroundControlPoint(row=4030.0, col=12660.0, x=99.47747019367358, y=-0.1356044843611649, z=-2.9550865292549133e-05, id='53', info=''), GroundControlPoint(row=4030.0, col=13926.0, x=99.36614759213732, y=-0.11163557152469691, z=-2.381112426519394e-05, id='54', info=''), GroundControlPoint(row=4030.0, col=15192.0, x=99.25483067731768, y=-0.08766771957079608, z=-1.829676330089569e-05, id='55', info=''), GroundControlPoint(row=4030.0, col=16458.0, x=99.14351900131565, y=-0.06370099042594742, z=-1.3002194464206696e-05, id='56', info=''), GroundControlPoint(row=4030.0, col=17724.0, x=99.03221217717538, y=-0.039735340073957244, z=-7.924623787403107e-06, id='57', info=''), GroundControlPoint(row=4030.0, col=18990.0, x=98.92090983767962, y=-0.01577079320831484, z=-3.07522714138031e-06, id='58', info=''), GroundControlPoint(row=4030.0, col=20256.0, x=98.80961163473663, y=0.00819257182000907, z=1.5581026673316956e-06, id='59', info=''), GroundControlPoint(row=4030.0, col=21522.0, x=98.69831726925788, y=0.03215478014871958, z=5.967915058135986e-06, id='60', info=''), GroundControlPoint(row=4030.0, col=22788.0, x=98.58702644036147, y=0.05611573340935314, z=1.0155141353607178e-05, id='61', info=''), GroundControlPoint(row=4030.0, col=24054.0, x=98.47573889490943, y=0.08007545677298648, z=1.4126300811767578e-05, id='62', info=''), GroundControlPoint(row=4030.0, col=25312.0, x=98.36515759039249, y=0.10388249215485533, z=1.784972846508026e-05, id='63', info=''), GroundControlPoint(row=6045.0, col=0.0, x=100.54289311574763, y=-0.5546486428260272, z=685.0006698500365, id='64', info=''), GroundControlPoint(row=6045.0, col=1266.0, x=100.43081958135876, y=-0.530502408840233, z=760.0007939795032, id='65', info=''), GroundControlPoint(row=6045.0, col=2532.0, x=100.32444503451912, y=-0.5075873813583428, z=425.00038781762123, id='66', info=''), GroundControlPoint(row=6045.0, col=3798.0, x=100.21289997027812, y=-0.4835519996737003, z=455.0004477193579, id='67', info=''), GroundControlPoint(row=6045.0, col=5064.0, x=100.10570348489681, y=-0.460456180774864, z=150.00005761999637, id='68', info=''), GroundControlPoint(row=6045.0, col=6330.0, x=99.99632394133126, y=-0.43688652482236146, z=-0.00013445783406496048, id='69', info=''), GroundControlPoint(row=6045.0, col=7596.0, x=99.8850117254185, y=-0.4128970800179826, z=-0.00012515205889940262, id='70', info=''), GroundControlPoint(row=6045.0, col=8862.0, x=99.77370592672115, y=-0.38890761152325504, z=-0.00011606886982917786, id='71', info=''), GroundControlPoint(row=6045.0, col=10128.0, x=99.66240603565707, y=-0.36491810864014773, z=-0.0001072026789188385, id='72', info=''), GroundControlPoint(row=6045.0, col=11394.0, x=99.55111157919816, y=-0.34092863301645115, z=-9.85562801361084e-05, id='73', info=''), GroundControlPoint(row=6045.0, col=12660.0, x=99.43982215357559, y=-0.31693913214645375, z=-9.013619273900986e-05, id='74', info=''), GroundControlPoint(row=6045.0, col=13926.0, x=99.32853736679876, y=-0.29294968505190516, z=-8.19377601146698e-05, id='75', info=''), GroundControlPoint(row=6045.0, col=15192.0, x=99.21725688710164, y=-0.2689602497401118, z=-7.396657019853592e-05, id='76', info=''), GroundControlPoint(row=6045.0, col=16458.0, x=99.1059803828889, y=-0.2449709280948571, z=-6.621330976486206e-05, id='77', info=''), GroundControlPoint(row=6045.0, col=17724.0, x=98.99470757843238, y=-0.22098168521175574, z=-5.8690086007118225e-05, id='78', info=''), GroundControlPoint(row=6045.0, col=18990.0, x=98.88343820487475, y=-0.19699256701941464, z=-5.1390379667282104e-05, id='79', info=''), GroundControlPoint(row=6045.0, col=20256.0, x=98.77217200171704, y=-0.17300368070970978, z=-4.4318847358226776e-05, id='80', info=''), GroundControlPoint(row=6045.0, col=21522.0, x=98.66090875270062, y=-0.1490150140907118, z=-3.7472695112228394e-05, id='81', info=''), GroundControlPoint(row=6045.0, col=22788.0, x=98.5496482317727, y=-0.12502668172399498, z=-3.0848197638988495e-05, id='82', info=''), GroundControlPoint(row=6045.0, col=24054.0, x=98.43839025341327, y=-0.10103867803568632, z=-2.445746213197708e-05, id='83', info=''), GroundControlPoint(row=6045.0, col=25312.0, x=98.32783765953045, y=-0.07720265243658453, z=-1.832377165555954e-05, id='84', info=''), GroundControlPoint(row=8060.0, col=0.0, x=100.50774341514176, y=-0.7367718994634289, z=480.00042406562716, id='85', info=''), GroundControlPoint(row=8060.0, col=1266.0, x=100.39983159662803, y=-0.7135056151048915, z=263.00012136437, id='86', info=''), GroundControlPoint(row=8060.0, col=2532.0, x=100.29109828634427, y=-0.6900591058200759, z=89.9998723231256, id='87', info=''), GroundControlPoint(row=8060.0, col=3798.0, x=100.1800971021884, y=-0.6661190664554266, z=74.99986377451569, id='88', info=''), GroundControlPoint(row=8060.0, col=5064.0, x=100.06986007049694, y=-0.642342246889066, z=-0.00024374108761548996, id='89', info=''), GroundControlPoint(row=8060.0, col=6330.0, x=99.95862392235972, y=-0.6183466281716569, z=-0.00023173633962869644, id='90', info=''), GroundControlPoint(row=8060.0, col=7596.0, x=99.84739029935825, y=-0.5943493069813199, z=-0.0002199467271566391, id='91', info=''), GroundControlPoint(row=8060.0, col=8862.0, x=99.7361590194588, y=-0.5703503975946956, z=-0.00020837318152189255, id='92', info=''), GroundControlPoint(row=8060.0, col=10128.0, x=99.62492993794879, y=-0.5463499081393922, z=-0.00019701942801475525, id='93', info=''), GroundControlPoint(row=8060.0, col=11394.0, x=99.51370289749377, y=-0.5223479681138787, z=-0.00018588360399007797, id='94', info=''), GroundControlPoint(row=8060.0, col=12660.0, x=99.40247776563571, y=-0.4983446428703032, z=-0.0001749703660607338, id='95', info=''), GroundControlPoint(row=8060.0, col=13926.0, x=99.291254442466, y=-0.4743398906688913, z=-0.00016428250819444656, id='96', info=''), GroundControlPoint(row=8060.0, col=15192.0, x=99.18003277422076, y=-0.4503339657730857, z=-0.00015381444245576859, id='97', info=''), GroundControlPoint(row=8060.0, col=16458.0, x=99.06881269265402, y=-0.42632675391283636, z=-0.00014357641339302063, id='98', info=''), GroundControlPoint(row=8060.0, col=17724.0, x=98.95759406878129, y=-0.40231846008786515, z=-0.000133562833070755, id='99', info=''), GroundControlPoint(row=8060.0, col=18990.0, x=98.846376806422, y=-0.3783091624157521, z=-0.00012377463281154633, id='100', info=''), GroundControlPoint(row=8060.0, col=20256.0, x=98.73516083959345, y=-0.3542988215931777, z=-0.00011421553790569305, id='101', info=''), GroundControlPoint(row=8060.0, col=21522.0, x=98.6239460447078, y=-0.33028769327163515, z=-0.00010489020496606827, id='102', info=''), GroundControlPoint(row=8060.0, col=22788.0, x=98.51273237751049, y=-0.30627567842088416, z=-9.578932076692581e-05, id='103', info=''), GroundControlPoint(row=8060.0, col=24054.0, x=98.40151973414879, y=-0.28226297742954837, z=-8.692033588886261e-05, id='104', info=''), GroundControlPoint(row=8060.0, col=25312.0, x=98.29101080409679, y=-0.25840137721294826, z=-7.833447307348251e-05, id='105', info=''), GroundControlPoint(row=10075.0, col=0.0, x=100.46939505009594, y=-0.9182065496957573, z=490.00042682141066, id='106', info=''), GroundControlPoint(row=10075.0, col=1266.0, x=100.36126640430022, y=-0.8948730063413666, z=290.00010570231825, id='107', info=''), GroundControlPoint(row=10075.0, col=2532.0, x=100.25418041798882, y=-0.8717625344140768, z=-0.0003975816071033478, id='108', info=''), GroundControlPoint(row=10075.0, col=3798.0, x=100.14295681782694, y=-0.8477519440844923, z=-0.0003826608881354332, id='109', info=''), GroundControlPoint(row=10075.0, col=5064.0, x=100.0317352883231, y=-0.8237387274977249, z=-0.00036793947219848633, id='110', info=''), GroundControlPoint(row=10075.0, col=6330.0, x=99.9205157454092, y=-0.7997228462570914, z=-0.00035342946648597717, id='111', info=''), GroundControlPoint(row=10075.0, col=7596.0, x=99.80929809085708, y=-0.7757043681993981, z=-0.0003391290083527565, id='112', info=''), GroundControlPoint(row=10075.0, col=8862.0, x=99.6980822171712, y=-0.751683438597612, z=-0.00032504182308912277, id='113', info=''), GroundControlPoint(row=10075.0, col=10128.0, x=99.58686803678636, y=-0.7276601374982283, z=-0.00031116604804992676, id='114', info=''), GroundControlPoint(row=10075.0, col=11394.0, x=99.47565548033354, y=-0.7036344843570027, z=-0.00029751006513834, id='115', info=''), GroundControlPoint(row=10075.0, col=12660.0, x=99.36444446990366, y=-0.6796065613351909, z=-0.0002840738743543625, id='116', info=''), GroundControlPoint(row=10075.0, col=13926.0, x=99.25323491846322, y=-0.6555765131282447, z=-0.00027085933834314346, id='117', info=''), GroundControlPoint(row=10075.0, col=15192.0, x=99.14202676950859, y=-0.6315443570440051, z=-0.0002578664571046829, id='118', info=''), GroundControlPoint(row=10075.0, col=16458.0, x=99.03081994042896, y=-0.6075102486417631, z=-0.0002451017498970032, id='119', info=''), GroundControlPoint(row=10075.0, col=17724.0, x=98.91961438138627, y=-0.5834742017342899, z=-0.00023255683481693268, id='120', info=''), GroundControlPoint(row=10075.0, col=18990.0, x=98.80841002884338, y=-0.5594363065634688, z=-0.00022024381905794144, id='121', info=''), GroundControlPoint(row=10075.0, col=20256.0, x=98.69720680933585, y=-0.5353967107200808, z=-0.00020816083997488022, id='122', info=''), GroundControlPoint(row=10075.0, col=21522.0, x=98.58600467742988, y=-0.5113554389303719, z=-0.00019630324095487595, id='123', info=''), GroundControlPoint(row=10075.0, col=22788.0, x=98.47400380528448, y=-0.4871385379320519, z=90.00002779252827, id='124', info=''), GroundControlPoint(row=10075.0, col=24054.0, x=98.36360343765513, y=-0.463268276420419, z=-0.0001732921227812767, id='125', info=''), GroundControlPoint(row=10075.0, col=25312.0, x=98.25310689563094, y=-0.4393745233880846, z=-0.00016220472753047943, id='126', info=''), GroundControlPoint(row=12090.0, col=0.0, x=100.4293937924664, y=-1.0992861993944323, z=610.0006292732432, id='127', info=''), GroundControlPoint(row=12090.0, col=1266.0, x=100.3271897112072, y=-1.0772176401926898, z=-0.0005689859390258789, id='128', info=''), GroundControlPoint(row=12090.0, col=2532.0, x=100.21597665248194, y=-1.05319117686908, z=-0.0005513597279787064, id='129', info=''), GroundControlPoint(row=12090.0, col=3798.0, x=100.10476519330317, y=-1.0291610748432545, z=-0.0005339309573173523, id='130', info=''), GroundControlPoint(row=12090.0, col=5064.0, x=99.99355530856856, y=-1.0051273730165213, z=-0.0005166986957192421, id='131', info=''), GroundControlPoint(row=12090.0, col=6330.0, x=99.88234694987946, y=-0.9810902195148092, z=-0.0004996638745069504, id='132', info=''), GroundControlPoint(row=12090.0, col=7596.0, x=99.77114008155343, y=-0.9570497006438244, z=-0.00048283208161592484, id='133', info=''), GroundControlPoint(row=12090.0, col=8862.0, x=99.65993467867264, y=-0.9330058504368318, z=-0.00046620890498161316, id='134', info=''), GroundControlPoint(row=12090.0, col=10128.0, x=99.54873070406501, y=-0.9089587598915964, z=-0.0004497962072491646, id='135', info=''), GroundControlPoint(row=12090.0, col=11394.0, x=99.43752810734121, y=-0.8849085816858656, z=-0.00043359212577342987, id='136', info=''), GroundControlPoint(row=12090.0, col=12660.0, x=99.32632686410072, y=-0.8608553450581583, z=-0.00041760317981243134, id='137', info=''), GroundControlPoint(row=12090.0, col=13926.0, x=99.21512692471181, y=-0.8367991978912817, z=-0.00040183402597904205, id='138', info=''), GroundControlPoint(row=12090.0, col=15192.0, x=99.10392826253225, y=-0.812740178903288, z=-0.0003862837329506874, id='139', info=''), GroundControlPoint(row=12090.0, col=16458.0, x=98.99273082868466, y=-0.7886784311796888, z=-0.00037095509469509125, id='140', info=''), GroundControlPoint(row=12090.0, col=17724.0, x=98.88153459628805, y=-0.7646139934228228, z=-0.0003558499738574028, id='141', info=''), GroundControlPoint(row=12090.0, col=18990.0, x=98.7703395293205, y=-0.7405469470332416, z=-0.00034097302705049515, id='142', info=''), GroundControlPoint(row=12090.0, col=20256.0, x=98.65914557560411, y=-0.716477449263897, z=-0.0003263242542743683, id='143', info=''), GroundControlPoint(row=12090.0, col=21522.0, x=98.54795271007686, y=-0.6924055293466905, z=-0.00031190458685159683, id='144', info=''), GroundControlPoint(row=12090.0, col=22788.0, x=98.436760881491, y=-0.6683313397700101, z=-0.00029771868139505386, id='145', info=''), GroundControlPoint(row=12090.0, col=24054.0, x=98.325570064647, y=-0.644254909774037, z=-0.0002837618812918663, id='146', info=''), GroundControlPoint(row=12090.0, col=25312.0, x=98.21508283803364, y=-0.6203284987521975, z=-0.0002701273187994957, id='147', info=''), GroundControlPoint(row=14105.0, col=0.0, x=100.40014349125441, y=-1.2827021840089496, z=-0.0007702391594648361, id='148', info=''), GroundControlPoint(row=14105.0, col=1266.0, x=100.28892609369387, y=-1.2586583886837548, z=-0.0007499316707253456, id='149', info=''), GroundControlPoint(row=14105.0, col=2532.0, x=100.17771063815958, y=-1.234610121845936, z=-0.0007298057898879051, id='150', info=''), GroundControlPoint(row=14105.0, col=3798.0, x=100.06649707543349, y=-1.2105575365067474, z=-0.000709855929017067, id='151', info=''), GroundControlPoint(row=14105.0, col=5064.0, x=99.95528536941288, y=-1.1865007239289107, z=-0.0006900951266288757, id='152', info=''), GroundControlPoint(row=14105.0, col=6330.0, x=99.84407549780676, y=-1.1624397088478713, z=-0.0006705224514007568, id='153', info=''), GroundControlPoint(row=14105.0, col=7596.0, x=99.73286742294721, y=-1.13837458720898, z=-0.0006511453539133072, id='154', info=''), GroundControlPoint(row=14105.0, col=8862.0, x=99.62166109482028, y=-1.1143055118810896, z=-0.0006319666281342506, id='155', info=''), GroundControlPoint(row=14105.0, col=10128.0, x=99.51045647625388, y=-1.090232573953198, z=-0.0006129862740635872, id='156', info=''), GroundControlPoint(row=14105.0, col=11394.0, x=99.39925355506132, y=-1.0661557458384388, z=-0.0005942154675722122, id='157', info=''), GroundControlPoint(row=14105.0, col=12660.0, x=99.2880522575846, y=-1.0420752894816396, z=-0.0005756514146924019, id='158', info=''), GroundControlPoint(row=14105.0, col=13926.0, x=99.17685256937143, y=-1.0179911867688414, z=-0.0005573006346821785, id='159', info=''), GroundControlPoint(row=14105.0, col=15192.0, x=99.06565445472464, y=-0.9939035192627953, z=-0.000539163127541542, id='160', info=''), GroundControlPoint(row=14105.0, col=16458.0, x=98.95445787378785, y=-0.9698123874888234, z=-0.0005212454125285149, id='161', info=''), GroundControlPoint(row=14105.0, col=17724.0, x=98.84326277751808, y=-0.9457179346314661, z=-0.0005035484209656715, id='162', info=''), GroundControlPoint(row=14105.0, col=18990.0, x=98.73206912683625, y=-0.9216202564154633, z=-0.0004860721528530121, id='163', info=''), GroundControlPoint(row=14105.0, col=20256.0, x=98.62087689669987, y=-0.8975193821923335, z=-0.0004688221961259842, id='164', info=''), GroundControlPoint(row=14105.0, col=21522.0, x=98.50968604894764, y=-0.8734154029776294, z=-0.0004517994821071625, id='165', info=''), GroundControlPoint(row=14105.0, col=22788.0, x=98.39849654538114, y=-0.84930840979027, z=-0.00043500587344169617, id='166', info=''), GroundControlPoint(row=14105.0, col=24054.0, x=98.28730833459231, y=-0.825198555227346, z=-0.00041844602674245834, id='167', info=''), GroundControlPoint(row=14105.0, col=25312.0, x=98.1768239878377, y=-0.8012382596044745, z=-0.00040221773087978363, id='168', info=''), GroundControlPoint(row=16120.0, col=0.0, x=100.36185107843961, y=-1.4641607670029828, z=-0.0009784381836652756, id='169', info=''), GroundControlPoint(row=16120.0, col=1266.0, x=100.25062990284555, y=-1.440095777274693, z=-0.0009556375443935394, id='170', info=''), GroundControlPoint(row=16120.0, col=2532.0, x=100.13941096406207, y=-1.4160256313903248, z=-0.0009329989552497864, id='171', info=''), GroundControlPoint(row=16120.0, col=3798.0, x=100.02819421296378, y=-1.3919504824892093, z=-0.0009105252102017403, id='172', info=''), GroundControlPoint(row=16120.0, col=5064.0, x=99.91697962759201, y=-1.3678703554634568, z=-0.0008882265537977219, id='173', info=''), GroundControlPoint(row=16120.0, col=6330.0, x=99.80576715656872, y=-1.3437854129213849, z=-0.0008661001920700073, id='174', info=''), GroundControlPoint(row=16120.0, col=7596.0, x=99.69455676524842, y=-1.3196957366514415, z=-0.0008441656827926636, id='175', info=''), GroundControlPoint(row=16120.0, col=8862.0, x=99.58334842875702, y=-1.2956013609248516, z=-0.0008224165067076683, id='176', info=''), GroundControlPoint(row=16120.0, col=10128.0, x=99.47214210993938, y=-1.2715023769791307, z=-0.0008008601143956184, id='177', info=''), GroundControlPoint(row=16120.0, col=11394.0, x=99.36093775839115, y=-1.247398937733212, z=-0.0007794974371790886, id='178', info=''), GroundControlPoint(row=16120.0, col=12660.0, x=99.24973534975024, y=-1.2232910726769795, z=-0.0007583387196063995, id='179', info=''), GroundControlPoint(row=16120.0, col=13926.0, x=99.13853483436095, y=-1.199178929942509, z=-0.0007373867556452751, id='180', info=''), GroundControlPoint(row=16120.0, col=15192.0, x=99.02733618560028, y=-1.1750625485054385, z=-0.0007166368886828423, id='181', info=''), GroundControlPoint(row=16120.0, col=16458.0, x=98.91523600130516, y=-1.1507449294128422, z=89.99960232060403, id='182', info=''), GroundControlPoint(row=16120.0, col=17724.0, x=98.80406346744186, y=-1.1266252464790398, z=89.99962767213583, id='183', info=''), GroundControlPoint(row=16120.0, col=18990.0, x=98.69303490192858, y=-1.1025326482483193, z=74.99960129894316, id='184', info=''), GroundControlPoint(row=16120.0, col=20256.0, x=98.58255944583387, y=-1.0785567054181757, z=-0.0006357943639159203, id='185', info=''), GroundControlPoint(row=16120.0, col=21522.0, x=98.47136954022021, y=-1.0544205925398324, z=-0.0006161350756883621, id='186', info=''), GroundControlPoint(row=16120.0, col=22788.0, x=98.36018126063391, y=-1.0302808446952307, z=-0.0005967002362012863, id='187', info=''), GroundControlPoint(row=16120.0, col=24054.0, x=98.24899459495502, y=-1.0061374298142371, z=-0.0005774945020675659, id='188', info=''), GroundControlPoint(row=16120.0, col=25312.0, x=98.13851206265464, y=-0.9821431713064229, z=-0.0005586380138993263, id='189', info=''), GroundControlPoint(row=18135.0, col=0.0, x=100.32352667375457, y=-1.645615836985871, z=-8.475035429000854e-08, id='190', info=''), GroundControlPoint(row=18135.0, col=1266.0, x=100.21230063400903, y=-1.621529577169001, z=-7.636845111846924e-08, id='191', info=''), GroundControlPoint(row=18135.0, col=2532.0, x=100.1010771275482, y=-1.5974374718335973, z=-7.35744833946228e-08, id='192', info=''), GroundControlPoint(row=18135.0, col=3798.0, x=99.98985610322565, y=-1.5733396837298688, z=-6.426125764846802e-08, id='193', info=''), GroundControlPoint(row=18135.0, col=5064.0, x=99.87863752604497, y=-1.5492362996125097, z=-6.426125764846802e-08, id='194', info=''), GroundControlPoint(row=18135.0, col=6330.0, x=99.76742137383519, y=-1.5251273444665, z=-5.8673322200775146e-08, id='195', info=''), GroundControlPoint(row=18135.0, col=7596.0, x=99.65620760900926, y=-1.5010129144862678, z=-5.4016709327697754e-08, id='196', info=''), GroundControlPoint(row=18135.0, col=8862.0, x=99.54499618058792, y=-1.4768931675367054, z=-5.029141902923584e-08, id='197', info=''), GroundControlPoint(row=18135.0, col=10128.0, x=99.43378705447641, y=-1.452768180718186, z=-4.7497451305389404e-08, id='198', info=''), GroundControlPoint(row=18135.0, col=11394.0, x=99.32258022661823, y=-1.428637943359777, z=-0.0009895442053675652, id='199', info=''), GroundControlPoint(row=18135.0, col=12660.0, x=99.21137560165732, y=-1.404502707966179, z=-0.0009657712653279305, id='200', info=''), GroundControlPoint(row=18135.0, col=13926.0, x=99.10017318997023, y=-1.380362397552697, z=-0.0009421929717063904, id='201', info=''), GroundControlPoint(row=18135.0, col=15192.0, x=98.98804584632798, y=-1.3560147715024389, z=89.99940622411668, id='202', info=''), GroundControlPoint(row=18135.0, col=16458.0, x=98.87687128514204, y=-1.331869839905853, z=89.99943534377962, id='203', info=''), GroundControlPoint(row=18135.0, col=17724.0, x=98.76569773921342, y=-1.3077199936848891, z=89.99946415517479, id='204', info=''), GroundControlPoint(row=18135.0, col=18990.0, x=98.65538469705602, y=-1.2837532802491698, z=-0.0008499100804328918, id='205', info=''), GroundControlPoint(row=18135.0, col=20256.0, x=98.54419269691113, y=-1.2595893722091387, z=-0.0008273618295788765, id='206', info=''), GroundControlPoint(row=18135.0, col=21522.0, x=98.43300264574611, y=-1.2354211075860828, z=-0.0008050361648201942, id='207', info=''), GroundControlPoint(row=18135.0, col=22788.0, x=98.32181454065453, y=-1.2112484117698752, z=-0.0007829293608665466, id='208', info=''), GroundControlPoint(row=18135.0, col=24054.0, x=98.21062832105052, y=-1.1870714803063211, z=-0.0007610432803630829, id='209', info=''), GroundControlPoint(row=18135.0, col=25312.0, x=98.10014653697382, y=-1.1630431847619476, z=-0.000739523209631443, id='210', info=''), GroundControlPoint(row=20150.0, col=0.0, x=100.28516976678081, y=-1.8270674267265, z=-1.126900315284729e-07, id='211', info=''), GroundControlPoint(row=20150.0, col=1266.0, x=100.17393778789366, y=-1.8029597585945378, z=-1.0523945093154907e-07, id='212', info=''), GroundControlPoint(row=20150.0, col=2532.0, x=100.06270861121182, y=-1.7788456886198296, z=-9.685754776000977e-08, id='213', info=''), GroundControlPoint(row=20150.0, col=3798.0, x=99.9514822269307, y=-1.754725184879666, z=-9.12696123123169e-08, id='214', info=''), GroundControlPoint(row=20150.0, col=5064.0, x=99.84025860015036, y=-1.7305983342880102, z=-8.568167686462402e-08, id='215', info=''), GroundControlPoint(row=20150.0, col=6330.0, x=99.72903767956794, y=-1.7064652997044918, z=-8.102506399154663e-08, id='216', info=''), GroundControlPoint(row=20150.0, col=7596.0, x=99.61781943064798, y=-1.6823261631780007, z=-7.450580596923828e-08, id='217', info=''), GroundControlPoint(row=20150.0, col=8862.0, x=99.50660384175761, y=-1.6581808975321024, z=-6.984919309616089e-08, id='218', info=''), GroundControlPoint(row=20150.0, col=10128.0, x=99.39539083651364, y=-1.6340297794188265, z=-6.51925802230835e-08, id='219', info=''), GroundControlPoint(row=20150.0, col=11394.0, x=99.2841804038943, y=-1.6098727768780086, z=-6.146728992462158e-08, id='220', info=''), GroundControlPoint(row=20150.0, col=12660.0, x=99.1721575582758, y=-1.5855318833181138, z=75.00028538610786, id='221', info=''), GroundControlPoint(row=20150.0, col=13926.0, x=99.06073053867254, y=-1.5613149234662476, z=98.00038053561002, id='222', info=''), GroundControlPoint(row=20150.0, col=15192.0, x=98.94963691031161, y=-1.5371646998407666, z=90.00035633612424, id='223', info=''), GroundControlPoint(row=20150.0, col=16458.0, x=98.83845995080367, y=-1.5129902740664412, z=90.00036310404539, id='224', info=''), GroundControlPoint(row=20150.0, col=17724.0, x=98.72816543999124, y=-1.4890030595400998, z=-4.470348358154297e-08, id='225', info=''), GroundControlPoint(row=20150.0, col=18990.0, x=98.61696962780452, y=-1.4648127895914338, z=-4.0046870708465576e-08, id='226', info=''), GroundControlPoint(row=20150.0, col=20256.0, x=98.50577611052873, y=-1.4406173369519806, z=-4.0046870708465576e-08, id='227', info=''), GroundControlPoint(row=20150.0, col=21522.0, x=98.3945848631235, y=-1.416416731423865, z=-3.5390257835388184e-08, id='228', info=''), GroundControlPoint(row=20150.0, col=22788.0, x=98.28339585619705, y=-1.3922110663980283, z=-0.0009938189759850502, id='229', info=''), GroundControlPoint(row=20150.0, col=24054.0, x=98.17220902013682, y=-1.3680004909374892, z=-0.00096922367811203, id='230', info=''), GroundControlPoint(row=20150.0, col=25312.0, x=98.06172691833218, y=-1.343938083651959, z=-0.0009450037032365799, id='231', info=''), GroundControlPoint(row=22165.0, col=0.0, x=100.24677984766224, y=-2.008515328147338, z=-1.4621764421463013e-07, id='232', info=''), GroundControlPoint(row=22165.0, col=1266.0, x=100.13554081614097, y=-1.98438629861146, z=-1.3783574104309082e-07, id='233', info=''), GroundControlPoint(row=22165.0, col=2532.0, x=100.024304909779, y=-1.960250059090999, z=-1.2852251529693604e-07, id='234', info=''), GroundControlPoint(row=22165.0, col=3798.0, x=99.91307208055436, y=-1.9361067583362, z=-1.210719347000122e-07, id='235', info=''), GroundControlPoint(row=22165.0, col=5064.0, x=99.80184229060949, y=-1.9119564976207994, z=-1.1455267667770386e-07, id='236', info=''), GroundControlPoint(row=22165.0, col=6330.0, x=99.69061551896357, y=-1.887799297452987, z=-1.0896474123001099e-07, id='237', info=''), GroundControlPoint(row=22165.0, col=7596.0, x=99.57939172715642, y=-1.863635259046927, z=-1.0151416063308716e-07, id='238', info=''), GroundControlPoint(row=22165.0, col=8862.0, x=99.46817086732183, y=-1.8394645262928193, z=-9.313225746154785e-08, id='239', info=''), GroundControlPoint(row=22165.0, col=10128.0, x=99.35695291455897, y=-1.8152871338511039, z=-8.940696716308594e-08, id='240', info=''), GroundControlPoint(row=22165.0, col=11394.0, x=99.2448999697859, y=-1.7909199957458697, z=75.00030375830829, id='241', info=''), GroundControlPoint(row=22165.0, col=12660.0, x=99.13354744466908, y=-1.76669891199205, z=90.00037243496627, id='242', info=''), GroundControlPoint(row=22165.0, col=13926.0, x=99.02227932059346, y=-1.7424892210340528, z=98.00041403621435, id='243', info=''), GroundControlPoint(row=22165.0, col=15192.0, x=98.91210933723987, y=-1.7185132265803247, z=-6.891787052154541e-08, id='244', info=''), GroundControlPoint(row=22165.0, col=16458.0, x=98.80090533377616, y=-1.6943040284774158, z=-6.51925802230835e-08, id='245', info=''), GroundControlPoint(row=22165.0, col=17724.0, x=98.6897039989066, y=-1.6700887819459758, z=-5.960464477539063e-08, id='246', info=''), GroundControlPoint(row=22165.0, col=18990.0, x=98.5785053067332, y=-1.6458675216475696, z=-5.774199962615967e-08, id='247', info=''), GroundControlPoint(row=22165.0, col=20256.0, x=98.46730920799399, y=-1.621640391305757, z=-5.4016709327697754e-08, id='248', info=''), GroundControlPoint(row=22165.0, col=21522.0, x=98.35611567664613, y=-1.59740742561451, z=-5.3085386753082275e-08, id='249', info=''), GroundControlPoint(row=22165.0, col=22788.0, x=98.2449246612732, y=-1.5731687777722207, z=-4.936009645462036e-08, id='250', info=''), GroundControlPoint(row=22165.0, col=24054.0, x=98.13373613767763, y=-1.5489244730174672, z=-4.563480615615845e-08, id='251', info=''), GroundControlPoint(row=22165.0, col=25312.0, x=98.02325265271054, y=-1.5248278787423222, z=-4.190951585769653e-08, id='252', info=''), GroundControlPoint(row=23071.0, col=0.0, x=100.22950777136563, y=-2.0900981424749836, z=-1.6298145055770874e-07, id='253', info=''), GroundControlPoint(row=23071.0, col=1266.0, x=100.11826337496667, y=-2.06595905490754, z=-1.5366822481155396e-07, id='254', info=''), GroundControlPoint(row=23071.0, col=2532.0, x=100.00702235836577, y=-2.041812475526032, z=-1.4621764421463013e-07, id='255', info=''), GroundControlPoint(row=23071.0, col=3798.0, x=99.89578466227567, y=-2.017658550718341, z=-1.3504177331924438e-07, id='256', info=''), GroundControlPoint(row=23071.0, col=5064.0, x=99.78455023892222, y=-1.99349737968575, z=-1.2759119272232056e-07, id='257', info=''), GroundControlPoint(row=23071.0, col=6330.0, x=99.67331905857438, y=-1.969328981114344, z=-1.1920928955078125e-07, id='258', info=''), GroundControlPoint(row=23071.0, col=7596.0, x=99.56209107501361, y=-1.945153454610789, z=-1.126900315284729e-07, id='259', info=''), GroundControlPoint(row=23071.0, col=8862.0, x=99.45086623347504, y=-1.9209709426437107, z=-1.0617077350616455e-07, id='260', info=''), GroundControlPoint(row=23071.0, col=10128.0, x=99.3396445029183, y=-1.8967814786147124, z=-9.96515154838562e-08, id='261', info=''), GroundControlPoint(row=23071.0, col=11394.0, x=99.2284258275192, y=-1.872585214727669, z=-9.592622518539429e-08, id='262', info=''), GroundControlPoint(row=23071.0, col=12660.0, x=99.11721016499457, y=-1.848382241576734, z=-9.033828973770142e-08, id='263', info=''), GroundControlPoint(row=23071.0, col=13926.0, x=99.00599750072791, y=-1.8241725217648606, z=-8.381903171539307e-08, id='264', info=''), GroundControlPoint(row=23071.0, col=15192.0, x=98.89478775161847, y=-1.7999563406101236, z=-7.916241884231567e-08, id='265', info=''), GroundControlPoint(row=23071.0, col=16458.0, x=98.78358091774507, y=-1.7757335944887653, z=-7.264316082000732e-08, id='266', info=''), GroundControlPoint(row=23071.0, col=17724.0, x=98.67237692956672, y=-1.7515045071486812, z=-6.984919309616089e-08, id='267', info=''), GroundControlPoint(row=23071.0, col=18990.0, x=98.56117574806913, y=-1.7272691600813423, z=-6.51925802230835e-08, id='268', info=''), GroundControlPoint(row=23071.0, col=20256.0, x=98.44997735771862, y=-1.7030275258269578, z=-6.332993507385254e-08, id='269', info=''), GroundControlPoint(row=23071.0, col=21522.0, x=98.33878169150391, y=-1.678779818763365, z=-5.774199962615967e-08, id='270', info=''), GroundControlPoint(row=23071.0, col=22788.0, x=98.22758870889655, y=-1.654526130040808, z=-5.494803190231323e-08, id='271', info=''), GroundControlPoint(row=23071.0, col=24054.0, x=98.1163983826215, y=-1.6302664892622194, z=-5.3085386753082275e-08, id='272', info=''), GroundControlPoint(row=23071.0, col=25312.0, x=98.00591327051733, y=-1.606154349054305, z=-5.029141902923584e-08, id='273', info='')], CRS.from_epsg(4326))

but yes, doing the interpolation properly with respect to the 3D terrain is another long discussion in itself.

The current stackstack behavior does seem confusing. 1. stackstac.stack could instead fail with "GCPs or lack of source CRS not supported". Or 2. since WarpedVRT is used by stackstack perhaps the basic GCP transform illustrated above could be used by default...

Learning towards solution 2. Would the fix then be to change this if-statement:

stackstac/stackstac/rio_reader.py

Lines 347 to 352 in 9106708

if self.spec.vrt_params != {
"crs": ds.crs.to_epsg(),
"transform": ds.transform,
"height": ds.height,
"width": ds.width,
}:

to allow for datasets with GCPs (but no CRSes)? Or maybe add an elif block?

@weiji14 weiji14 linked a pull request Sep 23, 2022 that will close this issue
Draft
@gjoseph92
Copy link
Owner

Yeah, agreed that overall stackstac could raise a better error here when ds.crs is None. (Note that no matter what, we don't know this until .compute() time when we actually fetch the data, so the error wouldn't come from stackstac.stack.)

I don't remember if warping from GCPs to a geotrans "just works" in GDAL, but I can confirm that with this one-line change, you no longer get an error:

diff --git a/stackstac/rio_reader.py b/stackstac/rio_reader.py
index 0e7c84a..b64427a 100644
--- a/stackstac/rio_reader.py
+++ b/stackstac/rio_reader.py
@@ -342,7 +342,7 @@ class AutoParallelRioReader:
                 )
 
             # Only make a VRT if the dataset doesn't match the spatial spec we want
-            if self.spec.vrt_params != {
+            if ds.crs is None or self.spec.vrt_params != {
                 "crs": ds.crs.to_epsg(),
                 "transform": ds.transform,
                 "height": ds.height,

Instead, you can see it "works":
Screen Shot 2022-09-23 at 12 16 51 PM

Of course, that data looks pretty meaningless. I'm not sure if that's because:

  1. that's what the GRD data actually looks like (certainly not true; see it in explorer)
  2. A WarpedVRT on a GCP doesn't "just work"
  3. stackstac is doing something else wrong (also think this is kinda unlikely)

It would be good to confirm this warping actually works before something like #182. Otherwise, it would be better to just raise an error than give back unusable data.

@weiji14
Copy link
Author

weiji14 commented Sep 23, 2022 •
edited
Loading

It would be good to confirm this warping actually works before something like #182. Otherwise, it would be better to just raise an error than give back unusable data.

Think we just need to pass in the correct epsg/crs parameters because the GCPs are in EPSG:4326 coordinates, but WarpedVRT doesn't know that, and it's still trying to reproject arbitrarily to EPSG:32647? Let me try around for a bit more.

Edit: Hmm yeah, tried setting src_crs in WarpedVRT at 4f656fe but still got the same purple/green/yellow output. The x/y coordinates look ok, I think it might be the resampling method that's the issue? Does anyone have another suitable example dataset (with gcps and no crs) that we can use to verify things?

@vincentsarago
Copy link

I'm not sure if this is still an issue but if you want to build a WarpedVRT with the gcps you need to do something like

import rasterio
from rasterio.transform import from_gcps
from rasterio.vrt import WarpedVRT

with rasterio.open("...") as src:
    with WarpedVRT(
        src,
        src_crs=src.gcps[1],
        src_transform=from_gcps(src.gcps[0]),
    ) as vrt:
        ...

https://github.com/cogeotiff/rio-tiler/blob/main/rio_tiler/io/rasterio.py#L102-L109

@weiji14
Copy link
Author

weiji14 commented Jan 6, 2023

I'm not sure if this is still an issue but if you want to build a WarpedVRT with the gcps you need to do something like

import rasterio
from rasterio.transform import from_gcps
from rasterio.vrt import WarpedVRT

with rasterio.open("...") as src:
    with WarpedVRT(
        src,
        src_crs=src.gcps[1],
        src_transform=from_gcps(src.gcps[0]),
    ) as vrt:
        ...

https://github.com/cogeotiff/rio-tiler/blob/main/rio_tiler/io/rasterio.py#L102-L109

Thanks @vincentsarago! I think we've tried setting src_crs, but not the src_transform yet. Haven't got the time to test this out now, but I've made a note at #182 (comment). Will try and pick this up next month again, unless someone wants to submit the bugfix sooner (happy for someone to take over #182).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Enable opening datasets with gcps even if no valid crs is present weiji14/stackstac
6 participants
@TomAugspurger @lossyrob @gjoseph92 @scottyhq @vincentsarago @weiji14