search

gjoseph92 / stackstac

Turn a STAC catalog into a dask-based xarray
https://stackstac.readthedocs.io
MIT License
238 stars 49 forks source link

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

Open weiji14 opened 1 year ago

weiji14 commented 1 year ago

Hi there, just encountering some issues while trying to stack Sentinel-1 GRD data from Planetary Computer at https://github.com/weiji14/zen3geo/pull/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:

```python-traceback --------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Input In [11], in () ----> 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 https://github.com/gjoseph92/stackstac/blob/9106708bcc20daec4c5975e9c1240de38c38f2f1/stackstac/rio_reader.py#L347-L362

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 :smile:

TomAugspurger commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 :slightly_smiling_face:

scottyhq commented 1 year ago

^ 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 https://github.com/corteva/rioxarray/issues/339#issuecomment-853877629 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 commented 1 year ago

@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:

https://github.com/gjoseph92/stackstac/blob/9106708bcc20daec4c5975e9c1240de38c38f2f1/stackstac/rio_reader.py#L347-L352

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

gjoseph92 commented 1 year ago

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 commented 1 year ago

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 4f656fe0c8fabcafa0ffeb301c7fb4f22a6af7cf 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 commented 1 year ago

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 commented 1 year ago

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 https://github.com/gjoseph92/stackstac/pull/182#discussion_r1063709499. Will try and pick this up next month again, unless someone wants to submit the bugfix sooner (happy for someone to take over #182).