1. Rasterio | Main
  2. Topics

Reprojecting GEOTIFF to epsg:3857 with scrambled CRS data (possibly GCJ-02)

toku2288@...
 

Hello, I am trying to reproject a population distribution raster file of China into 3857 space. I tried but the images are not aligning. Looking at the CRS it has a lot of different standards, I am not sure what to make of it.

The meta of the source file looks something like this:
{'driver': 'GTiff', 
'dtype': 'int32',
 'nodata': 2147483647.0,
 'width': 5049,
 'height': 5582,
 'count': 1,
 'crs': CRS.from_wkt('PROJCS["Albers_Conic_Equal_Area",GEOGCS["GCS_Unknown_datum_based_upon_the_Krassowsky_1940_ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krasovsky_1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'),
 'transform': Affine(1000.0, 0.0, -2638055.0614613006,
        0.0, -1000.0, 5922165.484718928)}
PopTIF.crs.data returns:
{'proj': 'aea',
 'lat_1': 25,
 'lat_2': 47,
 'lat_0': 0,
 'lon_0': 105,
 'x_0': 0,
 'y_0': 0,
 'ellps': 'krass',
 'units': 'm',
 'no_defs': True}
after running the reprojection algorithm:
{'driver': 'GTiff',
 'dtype': 'int32',
 'nodata': 2147483647.0,
 'width': 6708,
 'height': 6102,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:3857'),
 'transform': Affine(1184.7466003556895, 0.0, 7551104.353792087,
        0.0, -1184.7466003556895, 7299154.069454485)}




I am then using PopTIF.dataset_mask()  to mask the area of a Heightfield GEOTIFF (I requested the raster area of this file to be the same as the lat/lon of the PopTIF data:  'lat_1': 25,'lat_2': 47, 'lat_0': 0, 'lon_0': 105
 {'driver': 'GTiff', 
'dtype': 'float32',
 'nodata': -32767.0,
 'width': 2816,
 'height': 2288,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:4326'),
 'transform': Affine(0.02197265625, 0.0, 73.4765625,
        0.0, -0.02197265625, 53.7890625)}

 EPSG:4326 to EPSG:3857

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -32767.0,
 'width': 2594,
 'height': 2537,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:3857'),
 'transform': Affine(2655.775407243914, 0.0, 8179373.5227401415,
        0.0, -2655.775407243914, 7130308.041557059)}


So, After moving both GEOTIFFs into 3857 space. I then needed to upscale the Heightfield so the array matched the the size of the Population data:

import skimage.transform as st

pathToRaster  = r'./DEMTIF/cnl10_3857.tif'
newsize = reprojectedmask.shape


with rio.open(pathToRaster) as src:
    arr = src.read()
    arr[arr ==-9999] = np.nan
    arr  = arr[0,:,:]
    new_shape = newsize
    newarr= st.resize(arr, new_shape, mode='constant')

Resulting array can be masked, and looks something like this, Notice how there is no height data in areas such as Taiwan, There has been a mismatch with the projections, However I cant tell where/how. This may be due to the fact I specified the lat/lon in AEA(plus other CRS)? Would appreciate any guidance on how to find the correct area to sample and how to transform it 

toku2288@...
 

Hi, Just an update. Turns GeoTIFF does use GCJ encoding, However, There seem to be some open-source projects online which can reverse the transformation to high accuracy, Is there some way to use this algorithm to reproject the GeoTIFF file?

Even Rouault
 

On mardi 21 juillet 2020 18:52:10 CEST toku2288 via groups.io wrote:

> Hi, Just an update. Turns GeoTIFF does use GCJ encoding, However, There seem

> to be some open-source projects online which can reverse the transformation

> to high accuracy, Is there some way to use this algorithm to reproject the

> GeoTIFF file?

 

This had been discussed for PROJ in

https://github.com/OSGeo/PROJ/issues/982

 

Summary: apparently technically doable, but for legal reasons for folks operating in China, cannot be included in PROJ.

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com

1 - 3 of 3