Geoapplications development http://rgeo.wikience.org Higher School - PowerPoint PPT Presentation

Geoapplications development http://rgeo.wikience.org Higher School of Economics, Moscow, www.cs.hse.ru

2 Why do we use maps? Pictures from Лебедева О.А., Картографические проекции (методическое пособие), Новосибирский учебно - методический центр по ГИС и ДЗ, Новосибирск, 2000

3 Map projection (leads to distortions) • • • • • • • https://en.wikipedia.org/wiki/Geoid

4 Datum: aligning spheroid (ellipsoid) Geoid/Earth surface Spheroid 1 Spheroid 2

5 Putting all together

6 Notes on CRS definition select distinct coord_ref_sys_kind from epsg_coordinatereferencesystem

7 Notes on CRS definition (2) # R code # Convert (lat, lon) in meters to (lat, lon) in degrees library ( raster ) library ( rgdal ) # (lon, lat) in UTM projection – units are meters utm_c = cbind ( 470115, 6322515 ) utm_sp = SpatialPoints ( utm_c, proj4string = CRS ( "+init=epsg:32637" )) utm_sp # convert to (lon, lat) in WGS84 – units are degress global_sp <- spTransform ( utm_sp, CRS ( "+init=epsg:4326" )) global_sp

8 Ellipsoids • • •

9 Datum • • • https://en.wikipedia.org/wiki/World_Geodetic_System https://en.wikipedia.org/wiki/Geodetic_datum

10 Coordinate systems (CS) • • • •

11 Spherical geographic coordinate system Can you port that on Java from C++? http://www.codeproject.com/Arti cles/15659/Longitude-Latitude- String-Parser-and-Formatter ° ° ° °

12 Cartesian coordinate system

Projections Cylindrical Transverse Cylindrical Oblique Secant Cylindrical Cylindrical Conical Secant Conical Planar Secant Planar http%3A%2F%2Fwww.nps.gov%2Fgis%2Fgps%2F04_datums_coordinatesystems_65.ppt

14 Projections • • • • • • • http://www.nps.gov/gis/gps/04_datums_coordinatesystems_65.ppt

15 Why a map projection matters? 𝑦 1 − 𝑦 2 2 + 𝑧 1 − 𝑧 2 2 𝑒 =

16 Equirectangular projection Link 1

17 Equirectangular projection The plate carrée (French, for flat square ), is the special case where φ 1 is zero. • • • • • • Link 1

18 UTM: Universal Transverse Mercator • • http://www.geo.hunter.cuny.e du/~jochen/GTECH201/Lectur es/Lec6concepts/Map%20coor dinate%20systems/UTM%20a nd%20UPS.htm

19 UTM: Universal Transverse Mercator Apply a custom Transverse Mercator projection to each strip (zone) Pic 1 Pic2 Pic3

20 UTM: Universal Transverse Mercator • Custom means you have a separate cylinder for each strip ( zone ) • The cylinder is aligned to cross the spheroid along two lines ( in red ) • Each line is 180 km apart from the central meridian of a zone Pic 1 Pic2 Pic3 Pic4

21 UTM: zones • → • International Date Line) •

22 UTM: zones (another picture) • → • International Date Line) • PIC

23 UTM: N & S

24 UTM: 𝒛 False Easting and False Northing 𝒚 = 𝟒𝟏𝟏 𝟏𝟏𝟏 𝒛 = 𝟓 𝟏𝟏𝟏 𝟏𝟏𝟏 𝒚 = −𝟑𝟏𝟏 𝟏𝟏𝟏 without false easting False easting False northing North zones 500,000 m none South zones 500,000 m 10,000,000 m 𝒛 False northing – no negative y-coordinates for all locations in southern hemisphere zones (10,000,000 m is enough) False easting – no negative x-coordinates west of a zone's central meridian. Northern zones do not use false northing: 𝒚 their y-coordinates are naturally positive. Link 1 Link 2

25 UTM: distortion • • • Link 1

26 UTM: best practices • → • → Link 1

27 UTM: a limitation • • ° ° • Link 1

28 Universal Polar Stereographic system • • Link 1

29 Military Grid Reference System (MGRS) • • • •

30 Military Grid Reference System (MGRS) • Link 1 • o o o

31 NASA WW Example: MGRS Grid MGRS

32 Sinusoidal projection: very large scenes

33 Sinusoidal projection http://modis-atmos.gsfc.nasa.gov/MOD04_L2/grids.html https://en.wikipedia.org/wiki/Sinusoidal_projection

34 EPSG database: www.epsg.org Collection of CRS, CS, datum, etc. EPSG Dataset v8.7 contains: 5821 coordinate reference systems (CRS) from EPSG:200 to EPSG:69036405 126 coordinate systems (CS), CS != CRS 678 datum 14 variants of prime meridians 4232 CRS transformations (directly from → to) 1829 revisions ….. 21 tables

35 CRS Text Formats http://spatialreference.org/ GDAL (gdalsrsinfo): • PROJ.4 • OGC WKT • XML (GML based) • ESRI WKT format • Mapinfo style CoordSys format • simplified, etc. Try GeoTools, GDAL http://docs.geotools.org/stable/userguide/library/referencing/crs.html http://www.gdal.org/gdalsrsinfo.html

raro@ubuntu-pelligrini:/mnt/hgfs/RS_DATA/Landsat8/LC81790212015146-SC20150806075046$ 36 gdalinfo ./LC81790212015146LGN00_sr_band1.tif Driver: GTiff/GeoTIFF Files: ./LC81790212015146LGN00_sr_band1.tif GDAL output in WKT for Size is 8191, 8271 Coordinate System is: Landsat 8 Moscow scene PROJCS["WGS 84 / UTM zone 37N", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], don’t worry due to font size PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], – we are looking closer on it AUTHORITY["EPSG","4326"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], just in several slides PARAMETER["central_meridian",39], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","32637"]] Origin = (224385.000000000000000,6322515.000000000000000) Pixel Size = (30.000000000000000,-30.000000000000000) Metadata: AREA_OR_POINT=Area http://www.geoapi.org/3.0/javadoc/org/o Band_1=band 1 surface reflectance pengis/referencing/doc-files/WKT.html Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 224385.000, 6322515.000) ( 34d27'55.58"E, 56d57'49.93"N) Lower Left ( 224385.000, 6074385.000) ( 34d43' 3.37"E, 54d44'27.51"N) Upper Right ( 470115.000, 6322515.000) ( 38d30'26.82"E, 57d 2'42.39"N) Lower Right ( 470115.000, 6074385.000) ( 38d32' 5.80"E, 54d48'56.61"N) Center ( 347250.000, 6198450.000) ( 36d33'23.03"E, 55d54'25.94"N) Band 1 Block=8191x1 Type=Int16, ColorInterp=Gray Description = band 1 surface reflectance

37 EPSG codes interpretation One of the EPSG goals is to give short codes for frequent combinations to exchange metadata, the below WKT can be replaced by EPSG:32637 <projected cs> UTM 37 North PROJCS ["WGS 84 / UTM zone 37N", GEOGCS ["WGS 84", <geographic cs> DATUM ["WGS_1984", SPHEROID ["WGS 84", 6378137, 298.257223563, AUTHORITY[" EPSG "," 7030 "]], AUTHORITY[" EPSG "," 6326 "]], <prime meridian> PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY[" EPSG "," 4326 "]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",39], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT ["metre",1,AUTHORITY[" EPSG "," 9001 "]], AUTHORITY[" EPSG "," 32637 "]]

38 EPSG codes interpretation (2)

39 PROJ.4 gdalinfo -proj4 LC81790212015146LGN00_sr_band1.tif +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs The "no_defs" item ensures that no defaults are read from the defaults files. Sometimes they cause surprising problems. http://lists.osgeo.org/pipermail/mapserver-users/2003-November/046863.html

40 Notes on WGS84 As you can see from the EPSG database, WGS84 may mean (in different context) different entities: • Coordinate reference system • Datum • Ellipsoid

41 Ellipse (1) – ellipsoid (2) – spheroid (3) – oblate spheroid, ellipsoid of revolution, reference ellipsoid, Earth ellipsoid (4) • • • • Wikipedia: Spheroid Wikipedia: Ellipse • Wikipedia: Ellipsoid

42 Ellipse (1) – ellipsoid (2) – spheroid (3) – oblate spheroid, ellipsoid of revolution, reference ellipsoid, Earth ellipsoid (4) • ) • https://en.wikipedia.org/wiki/Earth_ellipsoid https://en.wikipedia.org/wiki/Reference_ellipsoid (2) (3) (4) SPHEROID ["WGS 84", 6378137, 298.257223563, AUTHORITY[" EPSG "," 7030 "]]

43 Reproject (raster=R, targetCRS=CRS, options) → →

44 Reprojection: the process http://demo.geo-solutions.it/share/OptimizingRasterReprojection.pdf

45 Reproject: visualization in web maps

46 Reproject use case: scenes from different projections