texts: GDAL upgrade causes errors in uploading shapefiles with ogr2ogr

texts:
After upgrading my GDAL version to 3.1.4, I am unable to upload shapefiles to a local postgres database by using ogr2ogr.  PostGIS is already installed and I am able to connect and run queries successfully on PGAdmin.  However, I am no longer able to upload shapefiles through the terminal.
The following command used to work:
ogr2ogr -f "PostgreSQL" PG:"host=localhost password=password dbname=postgres  port
=5432" USA_counties.shp -nln usa_counties

But now the following errors are thrown when I run it:
ERROR 1: ERROR:  could not load library "/usr/local/lib/postgresql/postgis-3.so": dlopen(/usr/local/lib/postgresql/postgis-3.so, 10): Symbol not found: _hash_bytes
  Referenced from: /usr/local/lib/postgresql/postgis-3.so
  Expected in: /usr/local/opt/postgresql/bin/postgres
 in /usr/local/lib/postgresql/postgis-3.so

ERROR 1: PROJ: proj_identify: Cannot find proj.db
ERROR 1: ERROR:  could not load library "/usr/local/lib/postgresql/postgis-3.so": dlopen(/usr/local/lib/postgresql/postgis-3.so, 10): Symbol not found: _hash_bytes
  Referenced from: /usr/local/lib/postgresql/postgis-3.so
  Expected in: /usr/local/opt/postgresql/bin/postgres
 in /usr/local/lib/postgresql/postgis-3.so

ERROR 1: ERROR:  could not load library "/usr/local/lib/postgresql/postgis-3.so": dlopen(/usr/local/lib/postgresql/postgis-3.so, 10): Symbol not found: _hash_bytes
  Referenced from: /usr/local/lib/postgresql/postgis-3.so
  Expected in: /usr/local/opt/postgresql/bin/postgres
 in /usr/local/lib/postgresql/postgis-3.so

ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db

Any ideas as to what steps need to be taken to successfully load postgis-3.so?
Edit:
When accessing
/usr/local/opt/postgresql/bin/postgres in the terminal, the following message appears:
postgres does not know where to find the server configuration file. You must specify the --config-file or -D invocation option or set the PGDATA environment variable.


texts: Exporting/Listing printed console results in a table in Google Earth Engine

texts:
I wrote a code with different function to bring me different output values e.g. the values at specific dates, statistics, cloud cover, etc. I printed them all individually in the console but now I would like to have one table or list with all the information so that I can export it easily.
Some of the values are in lists and some are in a feature collection.
I would like a table where for each date where there is an image within the given time, it shows the image name, the absolut value of point pt (B8 list), the standard deviation and the mean from statistics and the cloud cover score. The null values should also be included and just have a null values or empty cells.
It should conclude in a table/list with 36 dates.
How could I go about this?
var myB8 = l8.select("B8");

var getB8 = function(image) {

  // Reducing region and getting value
  var value_B8 = ee.Image(image)
    .reduceRegion(ee.Reducer.first(), pt)
    .get('B8');

  return value_B8;
};

var count = myB8.size();

var listOfImages = l8.toList(count);

var B8_list = listOfImages.map(getB8);

print("B8 list", B8_list);

var allDates = l8.aggregate_array('system:time_start');

var allDatesSimple = allDates.map(function procDates (ele) {

  return ee.Date(ele).format().slice(0,10);

});

var paired = allDatesSimple.zip(B8_list);

print (paired);

// RATIO
var RatioALL = idxList.map(function  calculateRatio (ele) {

  var idx = ele;

  var value = ee.Number(B8_list.get(idx));

  var image = listOfImages.get(idx);

  var ratio1 = ee.Image(image).select('B8')
                              .clip(poly)
                              .divide(value);

  var time = ee.Image(image).get('system:time_start');

  return ratio1.set('system:time_start', time);

});

print("Original Image Collection", myB8);

print ("Ratio All for Non Null Values", RatioALL);

//statistics
var getStats = function(image) {

  var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
  });

  var stats1 = ee.Image(image).reduceRegion({
    reducer: reducers,
    geometry: poly,
    scale: 15,
    bestEffort: true
  });

  return ee.Image(image).set(stats1);

};

var statistics = RatioALL.map(getStats);

print("statistics", statistics);

//Cloud cover
var getCloudScores = function(img){
    //Get the cloud cover
    var value = ee.Image(img).get('CLOUD_COVER');
    return ee.Feature(null, {'score': value})
};

var input = l8.map(function(image) { return image.clip(poly); });
var results = input.map(getCloudScores);
print ('Cloud cover', results);


texts: Listing active radio buttons within tab widget or layout?

texts:
I have a plugin widget with some tabs, and within the tabs a layout of radio buttons as shown below:

I would like to capture the list of active radio buttons within each tab like so:

Is there an easy way to do this dynamically  for each tab instead of programatically capturing the state of each radio button?

here is the structure of the UI in QT designer if it helps:


texts: How to export and import ST_GEOMETRY from Oracle environment to another Oracle environment?

texts:
I have a geographic layer stored in an Oracle Geodatabase with ST_GEOMETRY type.

I'd like to export that layer and its data into a SQL INSERT format in order to import it into another Oracle database which is in another server machine.

After the export I get lines like this:

Insert into myTable(OBJECTID,NAME,SHAPE) values ('2989','MyName','SDE.ST_GEOMETRY(8,1120,672587.528700002,6824864.9106,676521.2799,6829005.2801,NULL,NULL,NULL,NULL,8002536.39765853,16898.2791219904,3,'oracle.sql.BLOB@5e105008'));

I don't know how to deal with the BLOB type. I suppose oracle.sql.BLOB@5e105008references a memory address where some data is stored. How can I export that belob and import it in the new Oracle DB?


texts: Non-iterative algorithms for direct problem in the study of geodesic

texts:
I am finding an appropriate algorithm to implement in my program for direct problem in the study of geodesic. In brief, the direct problem deals with finding the end point (its longitude, latitude and azimuth) given the start point (its longitude, latitude and azimuth) and the distance from the start point to end point.

At the moment, I found out that Vincenty's formulae is a well-known algorithm for solving the geodesic direct problem. It is an iterative algorithm. According to this post, GeographicLib also has its own algorithms to solve the direct problem.

I have two questions. The first one is whether there are any non-iterative algorithms for this problem? If yes, can you give the name of the method or reference to it. The second one is if it is possible, can you give me the information about the algorithms used by GeographicLib.


texts: opengeo-suite upgraded from 3.0 to 3.1 on Ubuntu and I can't start GeoServer

texts:
I'm having trouble starting GeoServer after updating all of the opengeo-suite applications via apt-get update, apt-get upgrade and apt-get dist-upgrade.

GWC works, the dashboard is there as is geoexplorer. I can see geoserver in Tomcat's application manager, but I get a 404 tomcat error when I go to /geoserver.

Log file is here: http://pastebin.com/TdnYRxAZ

catalina.properties here: http://pastebin.com/acJCpAsa


texts: Convert multiple bsq files to jpg images in python

texts:
I am trying to convert a big number of bsq files to JPEG images. I found the way to do it only for one file but I cannot figure out how to put it in a loop in order to run repetitively all the files that are in a folder. Here is what I already did because I have already found it online.

import gdal

img_bsq  =image.bsq
img_jpg  =image.jpg

gdal.Translate(img_jpg, img_bsq, format='JPEG', width=1024, height=0,
scaleParams=[[]])


texts: What projection to use with Mapbox base/reference layers?

texts:
excuse me if this has been asked before, but I can not find answer after heavy googling.

What coordinate system/projection does shapefiles need to be in to be alligned correctly with mapbox reference layers?

This page states that the only projection supported is web mercator. However, when I project my shapefile to web mercator in ArcMap, the shapefile is displayed over Africa? Am I simply missing false easting and northing values?

Anyway, HERE are examples of alignments with two different projections.

As displayed above I've convert my shapefiles into different projections, without managing to align them with the reference layer. Help please?


texts: compile gdal with ESRI FileGDB_API Ubuntu

texts:
I'm trying to build GDAL on Ubuntu 16.04 with ESRI FileGDB_API driver support. I've followed the instructions here: https://trac.osgeo.org/gdal/wiki/FileGDB and have attempted to pass the location of proj included with the ESRI FileGDB_API. Here is my configure command:

./configure --with-static-proj=/usr/local/FileGDB_API/lib/libFileGDBAPI.so --with-fgdb=/usr/local/FileGDB_API

And the error:

configure: error: --without-static-proj and --with-fgdb are both available. There might be some incompatibility between system libproj and the embedded copy within libFileGDBAPI

I get the same error if I remove --with-static-proj and if I use --without-static-proj.

I'm attempting to compile gdal 2.3.1 and I'm using the FileGDB_API version 1.5 for gcc 5.1.


texts: When I get the total bounds of a shapefile, which are llcrnrlon, llcrnrlat, urcrnrlon, and urcrnrlat? Need it to use Maptlotlib Basemap

texts:
The goal is to use an ESRI webmap as a basemap beneath a shapefile for a script that visualizes earthquakes in the Philippines. But first, I'm practicing in using Basemap in a separate script.

According to the Basemap documentation, I need to get the following to dictate where the webmap is:

For most map projections, the map projection region can either be specified by setting these keywords:
Keyword     Description
llcrnrlon   longitude of lower left hand corner of the desired map domain (degrees).
llcrnrlat   latitude of lower left hand corner of the desired map domain (degrees).
urcrnrlon   longitude of upper right hand corner of the desired map domain (degrees).
urcrnrlat   latitude of upper right hand corner of the desired map domain (degrees).

So, to get the total bounds of the Philippines shapefile from GADM:

import geopandas as gpd
prov  = gpd.read_file('Provinces.shp')

In: prov.total_bounds
Out: array([116.9283371 ,   4.58693981, 126.60534668,  21.07014084])

To get an ESRI basemap:

import os
os.chdir(r'path')
import matplotlib.pyplot as plt
import numpy as np
os.environ['PROJ_LIB'] = r'path'
from mpl_toolkits.basemap import Basemap

tb = {'llon': 116.9283371, 'llat': 4.58693981,
      'ulon': 126.60534668, 'ulat': 21.07014084}

map = Basemap(projection="tmerc", ellps='WGS84', epsg='4326',
              llcrnrlon=tb['llon'], llcrnrlat=tb['llat'],
              urcrnrlon=tb['ulon'], urcrnrlat =tb['ulat'],
              lon_0=np.average([116.9283371, 126.60534668]),
              lat_0=np.average([4.58693981, 21.07014084]))

map.arcgisimage(service='ESRI_Imagery_World_2D', xpixels=2000,
                verbose=True, zorder=1)
plt.show()

Result is a blank plot

The code is inspired by this thread: Using Google Maps image as layer with Matplotlib Basemap?


texts: ArcGIS Linear referencing tools duplicating records and erroneous measures

texts:
I am using LocateFeaturesAlongRoutes_lr, MakeRouteEventLayer_lr, and MakeRouteEventLayer_lr to combine three different overlays into one based on this question I asked month ago. The code is as follow:

try:
    print "\n\t\t\tLocating features along roads",
    arcpy.LocateFeaturesAlongRoutes_lr(FTR_TO_RS_TEMP, ORIG_RD_FTR,
                                       "CONTROL_SE", "0 Meters",
                                       FTR_TO_RS_DB_TEMP,
                                       "RCSKEY LINE FMEAS TMEAS", "FIRST",
                                       "DISTANCE", "ZERO", "FIELDS",
                                       "M_DIRECTON")

    print "\t\t\t[OK]\n\t\tcombining features and PDT schedule",
    arcpy.OverlayRouteEvents_lr(FTR_TO_RS_DB_TEMP,
                                "RCSKEY LINE FMEAS TMEAS",
                                LP_SCHEDULE, "RCSKEY LINE START END",
                                "UNION", ROAD_LP_JOINED_TEMP,
                                "RCSKEY LINE FMEAS2 TMEAS2",
                                "NO_ZERO", "FIELDS", "INDEX")

    print "\t\t\t[OK]\n\t\t\tNew event route layer",
    arcpy.MakeRouteEventLayer_lr(ORIG_RD_FTR, "CONTROL_SE",
                                 ROAD_LP_JOINED_TEMP,
                                 "RCSKEY LINE FMEAS2 TMEAS2",
                                 LP_EVENT_LAYER_TEMP)

except Exception, e:
    tb = sys.exc_info()[2]
    print "Line %i" % tb.tb_lineno
    print e.message

For some reason, the routine is causing the following two issues:

Duplicate records of overlayed route events. For example, when ArcGIS creates a new record based on a union of two others, it would duplicate it 5 to 10 times.
New fmeas and tmeas that are -0.000023 to 0 and 0 to 2.999977 instead of 0 to 3 as in the original layers

I am quiet certain that the two issues are related and was wondering if anyone can shed a light


texts: How to convert different lidar data formats to raster DTM and measure its geographic position accuracy?

texts:
I'm trying to generate a DTM from lidar data and I'm new in the
this field.

I have different kind of files which describe the same area (different versions):

*.dtm
*.xyz
*.xyz with variable grid width (1 to 4 meter), according to the information sheet I have.

How would you handle each file format to get a DTM? Specially in the case of the XYZ-data with variable grid width, I have no idea about what to do.

Once you have the DTM, how can I be sure about the correct geographic position of the DTM?

I want to use the DTM to analyse the terrain. Which kind of format (of the DTM) would you recommend?

I use ArcGIS (with extensions) and would prefer a solution process on this platform. However, if there is an open source solution it's fine too.


texts: ArcGIS javascript API 3.6 How to prevent .on adding graphic layer from eliminating popups on other layers

texts:
I have am using mouse events to add a new graphic as such.

paCounties.on("mouse-over", function (evt) {
  var highlightGraphic = new Graphic(evt.graphic.geometry, highlightSymbol);
  map.graphics.add(highlightGraphic);
}

This works fine by showing my new graphic on mouse over.  However, I had three other layers that had popups and now the popups don't work.  If I comment out the above code and refresh my popups are back.  So something with adding the graphics layer or mouse-over event is preventing the popups from working on my other feature layers.

How can I have a mouse over and still retain popups on other feature layers?

EDIT: As per John's suggestion I tried using selection symbol but I can't seem to get it to work on mouseover.  I tried this:

I set my selection symbol on the feature layer.

paCounties.setSelectionSymbol(highlightSymbol);
paCounties.setRenderer(new SimpleRenderer(symbol));
map.addLayer(paCounties);

Then I try to trigger it on mouse over:

        paCounties.on("mouse-over", function (evt) {
            var query = new Query();
            query.geometry = evt.graphic.geometry;
            paCounties.selectFeatures(query, FeatureLayer.SELECTION_NEW);
        });

Which seems to highlight the county I am mousing over and all the counties around it.


texts: SVY and QDC files from Quickdraw Contours in QGIS

texts:
I am trying to view data collected from a Garmin sonar unit in QGIS.
The data on the Fishfinder shows or creates the contours but I can't seem to find a way to convert it to a format that QGIS can "display" or use.

The files are saved to the SD card as SVY file and as QDC file types. Does QGIS have anyway of even using these data sets?

Is there any software that can convert it into a usable file type that QGIS can use?


texts: Moving point feature and vertex from polyline at the same time

texts:
I have two layers in QGIS 3.14:

one with point geometry
one with polyline geometry

The point geometries (p) have the same coordinates as some of the polyline vertices (v):
(v)--------(v)----------(v/p)--------(v)-------------(v)----(v/p)-------(v/p)
Now I want to move a point geometry AND the vertex from the line at this point (v/p) to another coordinate:
(v)--------(v)-----------------(v/p)-(v)-------------(v)----(v/p)-------(v/p)
Is it possible to do this in one action? I think there must be a way to "connect" both layers, but I didn't find a function or a tool or a documentation for this.


texts: My RapidEye image won't appear in GRASS GIS

texts:
I have several RapidEye images that I want to use for object segmentation in GRASS GIS. When I import the the bands they all look like this and the same goes with a RGB composite:

I've used these images in ArcGIS and SAGA GIS before but never experienced this. The images are in tif. Any suggestions as to what the problem may be?


texts: How to read the saved extent of the QgsProject?

texts: I am trying to access the extent of QgsProject. How can i do that with PyQGIS?


texts: QGIS reshape tool and enable tracing (snap) tool

texts: I know how to increase or change the size of a polygon with the reshape tool, but can you also use Enable Tracing tool  with it, as you can with drawing a new feature?


texts: I have tiff file that has been modified and geotiff file with info, how do i convert my tiff to geotiff?

texts:
I have modified tiff file that i would like to turn to geotiff. My orginal geotiff has the following information that i got from gdalinfo:

Driver: Gtiff/Geotiff

Files: info.tiff

Size is 1127, 831 (new file is 4096 3020 in size)

Cordinate system is:

GEOGCS["WGS 84",

    DATUM[WGS_1984",

    SPHEROID[WGS 84", 6378137,298.257223563,

    AUTHORITY["EPSG","7030"]],

    AUTHORITY["EPGS","6326"]],

    PRIMEM["Greenwich",0],

    UNIT["degree",0,0174532925199433],

    AUTHORITY["EPSG","4326"]]

Orgin = (-17.385864000000002,28.226624000000001)

Pixel Size = (0.000277826086957, -0.000277632972323)

Metadata: AREA_OR_POINT=Area

Image Structure Metadata:

INTERLEAVE=BAND

Corner Coordinates:

Upper Left (-17.3858640, 28.2266240)

Lower Left (-17.3858640, 27.9959110)

Upper Right (-17.0727540, 28.2266240)

Lower Right (-17.0727540, 27.9959110)

Center (-17.2293090, 281112675

Band 1 Block=1127x3 type=int16, Colorinterp=gray

NoData Value=32768

Now i would like to place this new tiff that will be converted to same place with same information of the orginal data.

I think i need to use gdal_translate, i have seen some examples but those are too confusing to understand and i need some wrap thing too i think.

Would it be possible for someone to show me the lines i need to use with the provided information, so i could manage to do it by my self next time. I made the new image copy of the orginal, but with more resolution and slight tweaks, the new file is 4096 x 3020.

Thank You Very Much!


texts: Why are Landsat-8 derived NDVI values higher than SAVI values?

texts:
I have calculated NDVI and SAVI vegetation indices for a landsat-8 image. I am unexpectedly getting higher values for the NDVI than SAVI. What may be causing these unusual results?

NDVI=0.46
SAVI=0.34

Specific band pixel values are:

toa b4- 0.130
toa b5- 0.365
L= 0.5 taken


texts: Automating Model Builder (ArcMap) for NetCDF Processing

texts:
I have a NetCDF with 30+ years of monthly data from which I want to find the climatological normal. Using the NetCDF tutorial I found through NOAA (https://www.esrl.noaa.gov/psd/data/gridded_help/tools.html#arcgis in case anyone is interested), I discovered the Model Builder in ArcMap.

Is there anyway to have the model loop through each month's data to plug into the Raster Calculator so I can determine the climatological normal? I'd rather not have to build 360 branches manually. I can do this pretty easily in R, but I'm trying to expand my proficiency with Arc.


texts:
I use QGIS and I have a point and the road network. I need to automatically extract the road names that enclose the specific point.
Nearest neighbor analysis and buffer zones cannot do the job, because in many cases the point is closer, as a matter of measured distance, to neighbor roads and not the enclosing ones.
Are there any ideas on how only the enclosing roads could be extracted?


texts: Removing fields in an exported table using ArcPy

texts:
I am using ArcGIS Desktop 10.6: ArcMap 10.6. ArcPy script within a custom toolbox.
At the end of my script, I want to export the attribute table of a File Geodatabase Feature Class into a .csv file, however, I would like to only have specific fields within the attribute table to be exported into the .csv file.
I am getting the table to export, though, all the fields in the attribute table still remain and the one's meant to be removed are still there.
My only thought as to why this may not be working is because I am using a File Geodatabase Feature Class.


texts: Installing ArcGIS Python API with ArcGIS Pro

texts:
I am new to the world of ArcGIS and I come from a data science background. I recently installed ArcGIS Pro with a Single User license from my organization. I would like to use it with the Python API for deep learning.

Can anyone tell me how to get my setup properly working. I am confused as to how to login with the GIS() object in Python.

I would like to be able to perform tasks similar to this example: https://developers.arcgis.com/python/sample-notebooks/land-cover-classification-using-unet/


texts: How to get the geometry type of an empty PostGIS table?

texts:
I have a web application that needs to know the geometry type of a PostGIS table's geom field before inserting rows into it. I use the following query to determine the geometry type of the geom column:

SELECT GeometryType(geom) FROM my_schema.building LIMIT 1

This returns the geometry type of an actual row, so it does not work when my table is empty. How do I determine the geometry type of the geometry column itself?


texts: Convert shapefile to GeoJSON whilst stripping and renaming fields

texts:
I'm writing a Python script to convert shapefiles into GeoJSONs but I would like to remove many of the fields as they are redundant and rename the ones I am retaining.

Here's an extract of where I have got to so far.

ef convert_shape_to_geojson(indir, in_shape, in_projection):
    inShapefile = in_shape
    inDriver = ogr.GetDriverByName("ESRI Shapefile")
    inDataSource = inDriver.Open(inShapefile, 0)
    inLayer = inDataSource.GetLayer()

    outGeoSJON = 'output.geojson'
    outfile = os.path.join(indir, outGeoSJON)
    driver = ogr.GetDriverByName("GeoJSON")
    if os.path.exists(outfile):
        driver.DeleteDataSource(outfile)
    data_source = driver.CreateDataSource(outGeoSJON)

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(in_projection)

    out_layer = data_source.CreateLayer(outGeoSJON.split(".", 0), srs, geom_type=ogr.wkbMultiPolygon)

I've created new shapefiles in the past with brand new fields like so:

field_parcel = ogr.FieldDefn("PARCEL",ogr.OFTString)
field_parcel.SetWidth(20)
out_layer.CreateField(field_parcel)
out_layer.CreateField(ogr.FieldDefn("X",ogr.OFTInteger))
out_layer.CreateField(ogr.FieldDefn("Y",ogr.OFTInteger))
out_layer.CreateField(ogr.FieldDefn("BUFFDIS",ogr.OFTReal))

But I have no clue at the moment how to take existing field definitions from a source file but renaming the fields whilst also carrying over the source attributes.

Is anyone able to help out or point me to some further reading?


texts: How to zoom to specific lat/long point in QGIS

texts: Does QGIS have the ability to Zoom to a specific lat/long point upon the user entering one.  I can see it shows the coordinates in the bottom task bar but I'm not seeing a way to manually enter lat/long point and zoom to it.  Possibly even harder, I would like to zoom to specific lat/long point even though I'm in a UTM projection.


texts: qgis2web Leaflet - How to add the +/- symbol to reduce the layers in legend and how to add a logo in the webmap?

texts:
I'm newbie to edit HTML of Leaflet of qgis2web. I have some questions for you:

How can I put for each layer the sign "+" and "-" to reduce the legend of each layer.
How can I add a logo in the webmap?


texts:
I have just recieved 3.63TB of data for a project area and it is all in .sup and .smti file formats.

Looking up what these extensions are, I see it is a result of a product called SOCET SET (http://www.geospatialexploitationproducts.com/content/socet-set-v56/).

I am working in ArcGIS Desktop and see that support for .sup files has been revoked as of 10.2.

How do I begin to assess the data that I have in ArcGIS Desktop?


texts:
I published a table (File Geodatabase Table) to ArcGIS Online as a hosted feature for the purposes of creating an Operations Dashboard. Unfortunately, I have no idea how to update it.

I've been able to use the geoprocessing tool Replace Web Layer or the ArcGIS Pro Tool Replace Web Layer before for spatial datasets, but it specifies that my target layer (the hosted table) is unsupported or does not exist.

I've got a few workarounds I think could work, but I'd rather automate this process correctly by updating the hosted layer with my local data.


texts: Finding center of circle given three points?

texts:
I have been trying to solve this problem for quite some time, but can't find a good solution.

Given three sets of geographic coordinates how can I calculate the center point. By "center point" I mean the point on earth that is equidistant from the three points, on the same side of the earth. It would also be the circumcenter of a triangle formed by the three points.

This question is very similar to Calculate midpoint from a series of latitude and longitude coordinates, but that is dealing with a large number of points (of which all but three have been edited out) so the answers are either dealing with the centroid of a polygon (not what I'm looking for), are suggesting an average of the coordinates (not even close) or just don't answer at all. None of them are dealing with a circle.

I am dealing with circles that are no more than about 100 miles in radius. It doesn't need to be too ridiculously accurate, so an ellipsoid or even a spherical model would be plenty close enough. It would also work to transfer the coordinates to cartesian, solve the problem then transfer them back, but there is no easy way to do that given that the circle could be anywhere on earth. I'm not going to be working with areas anywhere near the poles, so that's not a consideration.

I am intending to use a spreadsheet to find this, but if there is some other resource that would be better I'm open to it.

Here are the things I'm finding online that are not helping:

People with the same problem keep being referred to the website geomidpoint.com. This site calculates your "personal center of gravity" from multiple points. I have no idea what it is calculating, but it is not the center of a circle with the points. In fact, more than three points are unlikely to be on the same circle.

There are plenty of math sites with solutions to this problem in cartesian coordinates, but that doesn't work with geographic coordinates.

There are a few places that are finding the geographic coordinates of the centroid of a triangle, but that is not the same as the circumcenter. The circumcenter of a triangle does not need to be inside the triangle.

There are a few questions over on different se sites that take a poke at it, but most math experts seem to have trouble applying it to geographic problems. That's why I'm posting here instead of on math.se.

Calculate the center point of multiple latitude/longitude coordinate pairs on SO
3D coordinates of circle center given three point on the circle on math.se

This guy on a math forum seems to have a good, but complex method. He's using a spherical model, so I don't know how far off it would be. I tried to follow his method, but it doesn't work out. I'm not sure if his formulas are wrong or if I'm doing something wrong in translating it to spreadsheet format.


texts: Not able to query on ArcGIS Server REST query end point

texts:
I have a map server published on ArcGIS Server and it has a Poles layer. It has following fields:

Now when I am trying to run any query against this layer by passing queries like

Owner IS NOT NULL

or

Owner = 'Utility'

I am getting this error while I am sure the Owner is not NULL and all has value of utility

  Failed to execute query.

Can you please let me know what I am doing wrong?


texts: How to show coincident point with different symbology in webgis?

texts:
our team have develop an webgis using ArcGIS silverlight API 3.1 .we have some points in database.These points are save in a table in sql server 2008R2 not in Geodatabase.The points have same longitudes and latitudes. I want to show these points with different symbols( colors or sizes) . I surf in the internet and found clustering . But the problem is that i have coincident points and clustering is use when the latitude and longitude of point are different. How can i show my points in webgis ?

Thanks for advance


texts: Merge several entities (polyline) into a single entity

texts:
Bonjour! j'ai une couche vecteur contenant des axes(itinéraire) routiers qui sont divisés chacun en des segments. j'aimerai fusionner ses segments en un seul polylignes par axe(itinéraire) et je cherche un outil qgis qui permet de le faire avec plusieurs segment.

I have a vector layer containing road axes (directions) which are each divided into segments. I would like to merge its segments into a single polylines per axis (route) and I am looking for a qgis tool which allows it to be done with several segments.


texts: Address and coordinates of European cities request

texts:
I have a database of couple of millions of records which I would like to geocode. The data consists of addresses and council. All addresses are of a certain European city.

Questions:

May I have accurate coordinates with a free/low cost geocoder? If yes, which one?
Can I buy from a company (like TomTom) for a fee, street address, number, postcode, council and accurate coordinates of European cities? Estimates of the fee price?


texts:
In what form should make requests to geoserver gml file from javascript?.

I have made this question, but neither on nor XMLRequest jquery and have downloaded some of geoserver, geoserver values always returns null.

someone will have examples, or know that library is useful for this case.

Thanks for your answers


texts:
Just as a thought idea.
I was thinking of creating a cell reception model using a cell towers point layer.
I want to create some sort of buffer radiating from each point but depending in other layers (buildings, surface roughness, etc..) I want the buffer to be longer or shorter in a specific direction.
Has anyone tried this?
Seen anything like this?
I am used to working with R and prefer creating such a tool with code and not a model in ArcGIS/QGIS if thats an issue.

Edit

so after being flagged as too broad, what i figure iam looking for is a way to create a buffer with a different width in every angle.

The point is being able to use some sort of formula (surface roughness, signal strength and so on..) to define the buffer's width in each angle.


texts: Make one of the Parameters an "Optional Parameter" in model - - ModelBuilder (ArcGIS Pro v2.5.0)

texts:
It is possible that there are several workarounds in ArcPy or scripts but I need this to work in Modelbuilder.

Simply I want the model to run with the populated and un-populated parameter fields (i.e. the Optional field is blank).

In my model I have the following parameters:
1. Street Address (point)(Required)
2. Property Boundary (polygon)(Required)
3. Buildings (polygon)(OPTIONAL).

The first step is to "Spatial Join" the Address to Properties (to apply the Address to the properties in my study area). Following this, I want to apply the address_property's address to the buildings that are within the property. However, sometimes I won't have the buildings feature and would like to continue with just the address_properties. Needless to say, this will form part of a larger tool.

The workaround to make this parameter Optional is to disconnect the input line and in the model's Properties -> parameters change it by the drop-down list. The thing is once all this is done and I try to run the model by right-clicking on the model and click on Open - the tool still asked me to populate the parameter.

This should be straight-forward but I possibly jinked it by saying that.


texts: PostGIS error reads "Geometry has Z dimension but column does not"

texts:
I try to store a polygon in GeoJSON format into my PostGIS table. Both are using the CRS EPSG:3857. This is how the polygon is defined. Note that I put in some whitespace to make it more readable.

{
    "type":"Polygon",
    "coordinates":[
        [-91.23046875,45.460130637921],
        [-79.8046875,49.837982453085],
        [-69.08203125,43.452918893555],
        [-88.2421875,32.694865977875],
        [-91.23046875,45.460130637921]
    ],
    "crs":{"type":"name","properties":{"name":"EPSG:3857"}}
}

To store that in my database I use the following SQL query.

INSERT INTO tablename (name, polygon)
VALUES ('name', ST_GeomFromGeoJSON('GeoJSON string comes here...'))

But for some reason the above query throws an error.

ERROR:  Geometry has Z dimension but column does not

What is wrong in my attempt?


texts:
I want to perform raster calculation from GDAL tool. I am using the formula as,

 gdal_calc.py -A Input_File.tif --outfile=yourOutput.tif -- calc="(A>=2200)*A" --NoDataValue=0 --type='Float32'

But the above statement gives an Error

  ImportError : numpy.core.multiarray failed to import.


texts: Transformation parameters for Everest datum to WGS84

texts:
I have a shapefile with the following projection information:

PROJCS["Nepal_india_boundary",GEOGCS["Nepal_india_Geo",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",500000.0],PARAMETER["Central_Meridian",84.0],PARAMETER["Standard_Parallel_1",26.666667],PARAMETER["Standard_Parallel_2",29.333333],PARAMETER["Scale_Factor",0.9997],PARAMETER["Latitude_Of_Origin",28.002808],UNIT["Meter",1.0]]

Which transformation parameters (preferably 7 parameters) can I use to transform the shapefile into WGS84 (decimal degree) system?


texts: How to differentiate areas which are spatially disconnected

texts:
So, I've done some geoprocessing on a bunch of raster data, and I have a large number of points spread out across a large area. I would like to remove all the data points which are only one or two cells, which is just noise, so that I can more easily identify areas where the data is continuous and relevant.

I'm thinking that the best way to do this is to convert to vector and separate each area which is not connected to another area into its own row, then delete all rows with a small area.

However, I don't know how to separate the data into many rows? Is there a way to do this?


texts: How do I json_encode without those quotes for the coordinates in leaflet geojson polygon?

texts:
I have a field name:  'coordinates'
I have the coordinates stored in there:  like this:
[[ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0] ]]

when I json_encode($geojson);

I get this output:

"coordinates": "[[ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0] "

notice the double quotes

it needs to be:

"coordinates": [[ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]

How do I json_encode without those quotes for the coordinates?

Here is my full code.

while ($row = mysql_fetch_assoc($rs)) {
    $properties = $row;
    $geometry = $row;
    # Remove fields
    unset($properties['coordinates']);
    unset($properties['type']);
    unset($geometry['unitID']);
    unset($geometry['milesSq']);
    unset($geometry['acres']);
    $feature = array(
        'type' => 'Feature',
        'properties' => $properties,
        'geometry' => array($geometry)
    );
    # Add feature arrays to feature collection array
    array_push($geojson['features'], $feature);
}

echo json_encode($geojson);

results:

{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"unitID":"001","milesSq":"111.000000","acres":"222.000000"},"geometry":[{"type":"Polygon","coordinates":"[[[-109.048595341028587, 40.834756771179649], [-108.966734749488324, 40.562222787215255],  [-109.048595341028587, 40.834756771179649]]]"}]}]}

ON JS:

$.ajax({
            type: "post",
            url: "<?php echo $this->app_lib->FOLDER['admin'].$this->router->class ?>/buscador",
            dataType: "json",
            data: {clat:center.lat, clng:center.lng, zoom:zoo, minlat:min.lat, minlng:min.lng, maxlat:max.lat, maxlng:max.lng, layers:<?php echo json_encode($layers); ?>},
            success: function(result){
              //if (data.error) {
                //map.removeLayer(data.points);
                $("#loader").hide()
                $.each(result, function(index, layer) {
                    if(layer.data.layer == 3){
                            $.each(layer.data.features, function(index, polygons) {
                                console.log(polygons.geometry.coordinates);
                                var poly = L.geoJson(campus, {
                                style: function (feature) {
                                    return feature.properties && feature.properties.style;
                                },
                                onEachFeature: onEachFeature,
                            });
                            map.addLayer(poly);
});


texts: How does ArcGIS calculate Live Traffic in a network dataset?

texts: I create a network analyst service in ArcGIS 10.2 with live traffic.  The traffic is classified into four classes (stop and go, slow, moderate and free flow). I know that when the service uses historical traffic the traffic calculate using free flow and travel speed.  In Live traffic I only have travel speed but the result of traffic in map shows the four classes.  How the service calculate and show live traffic in my service?


texts:
I am trying to update feature attributes from a .NET application that I'm working on.
I've cleaned the code down to the following example:

string serviceUrl = string.Format("http://***.cloudapp.net/arcgis/rest/services/{0}/{1}/FeatureServer", "***", "Test");
string featureUrl = string.Format("{0}/{1}", serviceUrl, "0");
var table = await ServiceFeatureTable.OpenAsync(new Uri(featureUrl));
if (table.IsInitialized)
{
var found = await table.QueryAsync(1);
    if (found != null)
    {
        found.Attributes["testtext"] = "test";
        await table.UpdateAsync(found);
        if (table.HasEdits)
        {
            var result = await table.ApplyEditsAsync();
        }
    }
}

It finds the feature with OBJECTID = 1, table.HasEdits is true after the update, and the result variable contains a ResultItem in the UpdateResults property, with no errors and Success = true.

So no errors and no indication something went wrong, but the attribute table in ArcGIS doesn't change at all. If I try to change the geometry of the feature, by adding a line like this:

found.Geometry = new Polygon(new List<MapPoint> { new MapPoint(0, 0), new MapPoint(-5, 0), new MapPoint(-5, -5), new MapPoint(0, -5), new MapPoint(0, 0) });

It will update the geometry correctly... but still no update to the other attributes of the feature.

The feature class in question here, has a field called "testtext" which is of type "Text". I've also tried fields with type short, double, date etc, but I'm still not able to update or add any attribute data.

Any clues as to what I'm doing wrong here?


texts: Filling no data in raster using focal statistics in batch in ArcMap

texts:
I know how to use the below python script to fill no data in a single raster using focal statistics, but I have a series of rasters, I do not want to run the script one by one.
Can any people tell me how to do it in batch?
>>>
... import arcpy

... from arcpy import env

... from arcpy.sa import *

... env.workspace = "C:/Users/wing/Desktop/try/NDVI"

... outCon = Con(Raster("MOD13Q1.A2000049.250m_16_days_NDVI.tif")<-2000, FocalStatistics("MOD13Q1.A2000049.250m_16_days_NDVI.tif",NbrRectangle(3,3),"MEAN","NODATA"), "MOD13Q1.A2000049.250m_16_days_NDVI.tif")

... outCon.save("C:/Users/wing/Desktop/try/outcon")


texts: How to calculate elevation of the point?

texts:
I have two points: A ( x1 ( latitude ), y1 ( longitude ), z1 ( elevation ) ) and B ( x2, y2, z2 ). x1,y1,z1 and x2,y2 are known, also I know the distance between A and B.
How can I calculate z2 without using internet, totally offline ? ( I am using Google Maps on Android )


texts: Implementing transactions with Collector for ArcGIS?

texts: Currently I am looking into using ArcGIS Collector for field use.  I would like for there to be a QA process done before final submission. Is there a way to implement an approval process for the submissions similar to SQL transactions?


texts: Getting Coordinates from a feature with VBA

texts:
Is there a way to get the longitude and latitude coordinates from a feature class? I have some code which can select the first layer in the Table of Contents. I want to be able to get the lat./long. of a feature in that layer (any lat./long. within the Data Frame's current extent will do).

With a point file I can see this being pretty straightforward (using the IPoint interface), but how would I go about do this with a polygon?

(Bonus points to anyone who can tell me if I can get the lat./long. using a FactoryCode!)

Here's my code thus far:

Dim pMxDoc As IMxDocument
Dim pMap As IMap
Dim pTargetFL As IFeatureLayer
Dim pTargetGD As IGeoDataset
Set pMxDoc = ThisDocument
Set pMap = pMxDoc.FocusMap
Set pTargetFL = pMap.Layer(0) 'Set the IFeatureLayer as the first layer in the TOC
Set pTargetGD = pTargetFL.FeatureClass

Dim pTargetSR As ISpatialReference
Dim pGeoDatasetEdit As IGeoDatasetSchemaEdit
Dim pSpatRefFact As ISpatialReferenceFactory
Dim pProjCoordSys As IProjectedCoordinateSystem
Set pTargetSR = pTargetGD.SpatialReference 'Set the spatial reference object

Dim fatCode As Integer
fatCode = pTargetSR.FactoryCode
MsgBox ("Layer Name: " & pTargetFL.Name & vbNewLine & "Layer Projection: " & pTargetSR.Name & vbNewLine & "Factory Code: " & fatCode)


texts: Enabling Surface Constraints on a LASD in ModelBuilder?

texts:
I'm working on a LAS processing model. I have been able to:

create a LASD from a folder of LAS files.
create a mask from the LASD.
apply the mask (called Surface Constraints) to the LASD file.

However, because I cannot toggle it as "selected" via a ModelBuilder tool, then none of the subsequent analysis honors the mask. I have checked that running my next process (LASD to Raster) does behave as expected if I select it.

Am I missing a tool or is it missing from the tool set for 10.3? If it's available in any of the upgrades I'd be interested to know (10.4, 10.5, Pro 1.4) as well.


texts:
Since QGIS 2.0.1 the default color ramps and styling symbols seem to disappear for some users.
For me, this started with QGIS 2.0.1 and still consists with 2.4.
I can neither choose Smybols when styling a line vector layer, nor choose from default color ramps.

I'm looking for a universal solution to tackle this problem. There is already an open bug report, but it seems to be forgotten.


texts: Installing packages for ACOLITE atmospheric correction; RuntimeError - HDF4 DLLs

texts:
I would like to use ACOLITE atmospheric correction with python, however, I cannot install pyhdf.
The RuntimeError is following:
RuntimeError: Cannot find required HDF4 DLLs -- check LIBRARY_DIRS

I was googling, however I didn't find any specific answer for the error.
What should I do?
I am a Windows user.  However, I can use in an emergency the WSL Debian.


texts: Format Issue with Dynamic Text using Data-Driven Pages

texts:
I have a dynamic text element I've edited to show open and closed store locations:

I cannot seem to remove the spaces between my text. Note: the line you see before the last "0" is my legend border. The proper display should look like (not bold):
2014: 3/1 2015: 0/1 2016: 0/0

Here is my current text properties:

2014:<dyn type="page" property="attribute" field="StoreB_Cls2014" domainlookup="true"/>/<dyn type="page" property="attribute" field="StoreA_Cls2014" domainlookup="true"/>2015:<dyn type="page" property="attribute" field="StoreB_Cls2015" domainlookup="true"/>/<dyn type="page" property="attribute" field="StoreA_Cls2015" domainlookup="true"/>2016:<dyn type="page" property="attribute" field="StoreB_Cls2016" domainlookup="true"/>/<dyn type="page" property="attribute" field="StoreA_Cls2016" domainlookup="true"/>

The attributes are Double format.


texts: GeoTools modify a FeatureType to change its name

texts:
I have a featureCollection and I want to save it into PostGIS dataStore. So I'm going to create a database table by geotools. I want to name my table "ali". But I can't rename the schema!

DataStore dataStore = DataStoreFinder.getDataStore(params);
SimpleFeatureCollection collection=...;
SimpleFeatureType schema = collection.getSchema();
SimpleFeatureType schema2 = DataUtilities.createSubType(schema,null,null,"ali",null);
dataStore.createSchema(schema2); // but it creates with name "polygon"

How can I rename the schema to create a table with my desired name?


texts: Modifying “pg_hba.conf” file to listen a dynamic address

texts:
this is the sequel of "Modifying "pg_hba.conf" file to listen local TCP/IP conection"

I manage to modify pg_hba.conf to listen other PCs and allow them through the firewall, but my problem now is that the addresses are dynamic, so the configuration work for one day, but in the next day the router assigns the addresses randomly and we are in trouble again. When I use mySQL I use to work it around using the machine name as the server, but if I use it in the pg_hba it doesn't work. I show you:

what works:

# TYPE  DATABASE        USER            ADDRESS                 METHOD
# IPv4 local connections:
  host    all             all             192.168.1.18/32         md5

what doesn't:

# TYPE  DATABASE        USER            ADDRESS                 METHOD
# IPv4 local connections:
  host    all             all             the_mch_name/32         md5

Is there any solution besides static IP?


texts: How to use an ArcGIS Server globe service on the web?

texts:
I created an ArcGlobe document containing a 3D multipatch building...

... then published this to ArcGIS Server as a Globe Service. The help page Displaying data from a globe service says:

  Globe services can be consumed and displayed by the suite of ArcGIS
  applications, including ArcGlobe, ArcCatalog, ArcReader, and ArcGIS
  Explorer Desktop

but I can't see any mention of a web API or platform. What is the next step in order to use this globe service in a web page?


texts: Calculate histogram from raster stored in GeoServer

texts:
I have a raster stored in GeoServer which I can show it with Openlayers in my web page:

raster = new TileLayer({
  source: new TileWMS({
    url: 'http://myip/geoserver/myworkspace/wms',
    params: {'LAYERS': 'myworkspace:mylayer','TILED': true},
    serverType: 'geoserver',
  });
});

But I want to show also an histogram. I found a way to get the value of a pixel given a coordinate with the following url:

urlRequest = this.raster.getSource().getGetFeatureInfoUrl(
    myCoordinate,myview.getResolution(),
    'EPSG:3857',{'INFO_FORMAT': 'application/json'});

Is there a similar request for getting all values from the raster given an extent instead of a single coordinate or is there a better way to calculate an histogram?


texts: Apply Database based Rules to Openlayers

texts:
im trying to apply rules, where the data is fetched from an external database.
This is my code:

// create a new style and add rules
            var defStyle = new OpenLayers.Style({
                strokeWidth: 1
            });

        defStyle.addRules([

            new OpenLayers.Rule({
                symbolizer: {
                    fillColor: "#ffeecc"
                },
                filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.LESS_THAN,
                    property: "order07",
                    value: "04000"
                })
            }),

            new OpenLayers.Rule({
                symbolizer: {
                    fillColor: "#ffcc99"
                },
                filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.BETWEEN,
                    property: "order07",
                    lowerBoundary: "04000",
                    upperBoundary: "08000"
                })
            }),
            new OpenLayers.Rule({
                symbolizer: {
                    fillColor: "#ff9966"
                },
                filter: new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.BETWEEN,
                    property: "order07",
                    lowerBoundary: "08000",
                    upperBoundary: "11000"
                })
            }),
        ]);

        var selStyle = new OpenLayers.Style({
            fillColor: "#ffaa00"
        });

        // combine styles in a style map
        var myStyleMap = new OpenLayers.StyleMap({
            "default": defStyle,
            "select": selStyle
        });

what i want to do is something like this:

I have a database with

ID| Name

1 | Bob

2 | Sam

and i want the rules getting applied according to one of those values.
Is it possible to include external database information that is not in the shapefile  ??

In spoken words the original rule is "Color all vectors where the property:order07 is less than value:04000 with the color #ffeecc"

and i want to change it into something like "Color the vectors where the property:fromDatabase xy: is value:xxyy with the color #xxxxxx"


texts: Linking Access Database with arcgis 10 and is it more efficient?

texts:
I have currently started working for a fisheries consultancy. They have all their data eg. catch rates, fishing effort, gear types in an Access database. Existing and previous projects have map/mxd files saved on the c drive.

Requests to edit  an existing map will come through e.g update this map with new fishing effort data. They currently run a query in Access, then import the result into excel and create a table, they then finally add the table to arcgis desktop 10 and then convert into a shapefile by extracting xy coordinates.

Wouldn't it be simpler to link the Access database to arc and select the data required for the update directly making this mundain task quicker?

And isn't it simple to run a access like data query within arc?


texts: How to get coordinates of Drag Box using OpenLayers 3?

texts:
I'm pretty new to From Openlayers and I have a problem. I can draw a rectangle using DragBox on my OpenLayers 3 map. How can I get its coordinates? I tried getGeometry() but it returns 'null'. I think I just added the 'drawing dragbox feature' to my map and I am not creating an actual rectangle. I am going to use the area inside these coordinates later. Here is my code:

        boxControl =new ol.interaction.DragBox
        ({
            condition: ol.events.condition.altKeyOnly,
            style: new ol.style.Style
            ({
                stroke: new ol.style.Stroke({color: [0, 0, 255, 1]})
            })
        });
        map.addInteraction(boxControl);

OpenLayers 3 is pretty new and most of the resources on the web are for 2.x and I couldn't find any help.


texts: OpenGeo SDK modify the getfeatureinfo popup

texts:
I would like to modify the opengeo sdk gxp.FeatureEditPopup.
I want to construct a clickable hyperlink, I would like to add a new field to the popup that would be "http://someurl.com/" followed by a field from my data: "http://someurl.com/someDataField".
Any ideas would be appreciated.
Ben


texts: Postgis, Installing Tiger Geocoder :Using Extension

texts:
I am trying to install the tiger geocoder for postGIS. I am using the 2.2 version of postGIS on Windows 10. I am using this set of instruction to do so.

I do not know the meaning of the fourth step.

I am confused by what they mean be edit the paths in the declare_sect column to those that fit your pg, unzip,shp2pgsql, psql, etc path locations.

Here's the code that I am using to do this.

INSERT INTO tiger.loader_platform(os, declare_sect, pgbin, wget, unzip_command, psql, path_sep,
           loader, environ_set_command, county_process_command)
SELECT 'debbie', declare_sect, pgbin, wget, unzip_command, psql, path_sep,
       loader, environ_set_command, county_process_command
  FROM tiger.loader_platform
  WHERE os = 'sh';

Can anybody tell me what to do. Please write down the specific steps.


texts:
As in the Hillshade Example, I wanted to automate adding all Layer from one Image for visual inspection.

var band_names =Image.bandNames()
print(ee.String(band_names.get(0)))
  for (var i = 0; i < 4; i += 1) {
  Map.addLayer(Image.select(ee.String(band_names.get(i))), {}, ee.String(band_names.get(i)));
}

This returns the longest Error ever, not gonna post it entirely

2 Questions:

the Layer name is a String:

print(ee.String(band_names.get(0)))

yields a string. What's wrong with this code?
How do I call the number of Bands for the automation? I imagine something like :

ee.Number(ee.Image.bad


texts: Finding coordinates of point on line joining two points at certain distance from one of points using shapely

texts:
I am new to GIS. The question I am going to ask has been asked before but it didn't solve my purpose. I am only using python and shapely. I need to find the coordinates of a point on a line joining two points at a given distance from one of the points. To help illustrate, I am including the following figure:

In the figure, I need to find the coordinates (lat/lon) values of the red points at a certain distance (say d1*1.01) from the point marked as C. Given are the points C, and P1,P2 and P3 with their coordinates in lat/lon format. How do I do find these values (lat/lon values of the red points) using only python and shapely. I don't want to incorporate Arcpy unless there is no other options.


texts: Dissolve fails, geometry errors, and cleaning problems

texts:
I fairly new to GIS, and I'm having an issue with dissolving a vector. To start, I'm running QGIS 2.14.15 on Windows 7 Pro. When I try to dissolve the polygons it creates a blank shapefile. All of the files are in the same projections. I poked around and found Fixing geometry validity errors in QGIS? which lead me to finding two errors (see below).

187,Unknown geometry type 0;
187,Geometry has 1 errors.
221,Unknown geometry type 0;
221,Geometry has 1 errors.

When I select the an entry and try to Pan to selection it says "Geometry is Null". Running LWGEOM - Make Valid gives "[Error 87] The parameter is incorrect See log for more details", but there aren't any details in the log. Grass 7 v.clean runs and doesn't seem to do anything (not sure if I'm using it correctly).

Edit: The files are set to EPSG:3310. I used union to join two polygon layers showing range distribution. My dissolve process argument is ---Dissolve all---.  When using the check geometry validity function on the union layer, I selected Save Error Location as a point shapefile, the file is empty. In LWGEOM I used the default setting. Using Grass v.clean I used cleaning tools snap, bpol, prune and rmarea. I systematically varied the threshold between 0.1 and 1000. In the advanced parameters varied v.in.ogr snap tolerance from -1.0 to 5. and the min area from 0.001 to 1000. They all give the same empty file.

Since my OP I tried using SAGA 2.3.2 => Vector Polygon Tools=> Polygon dissolve (all polygons). It will dissolve the layer but omits 5 polygons where data intersects.


texts: Call WFS request after the map load in OpenLayers 3

texts:
I am trying to call WFS service and store the response in a variable after the map is loaded completely in OpenLayers 3. Below is my code

this.map = new ol.Map({
        interactions: olgm.interaction.defaults(),
        controls: ol.control.defaults({
            attribution: false,
            zoom:false
        }).extend([
            new ol.control.ScaleLine()
        ]),
        target: olMapDiv,
        view: view
    });

this.map.once('postrender', function(event){
        self.test();

    });

this.test = function () {
        self.layerWFS = new ol.layer.Vector({
          source: new ol.source.Vector({
            loader: function(extent) {
              $.ajax('http://localhost:8080/geoserver/geowebmap_ws/ows', {
                type: 'GET',
                data: {
                  service: 'WFS',
                  version: '1.1.0',
                  request: 'GetFeature',
                  typename:'geowebmap_ws:tn_highways',
                  srsname: 'EPSG:3857',
                  outputFormat: 'application/json',
                  bbox: extent.join(',') + ',EPSG:3857'
                }
              }).done(function(response) {
                    self.layerWFS
                    .getSource()
                    .addFeatures(new ol.format.WFS()
                      .readFeatures(response));
                  });
            },
            strategy: ol.loadingstrategy.bbox,
            projection: 'EPSG:3857'
          })
        });
            mapService.layer_wfs = self.layerWFS;
            console.log(self.layerWFS);
    }

Am not getting any WFS features in self.layerWFS after the map is loaded. Am I missing anything?


texts: Identify false outlets on a river network

texts:
I have a river network that is very poorly defined (lots of gaps).  I need to identify all the false outlets.  The picture below provides an example of one of these false outlets.  The image shows that the flow path (yellow) ends ubrubtly (circled in red), it is these end point I need to identify.

I can already identify all the outlets (last point in a reach that is not a junction) but I'm not sure how I approach defining if it is a true outlet or a false outlet. I have access to python, QGIS, GRASS and SAGA.


texts: Is it possible to stop an algorithm midway through processing in QGIS?

texts: This question come about do to an annoying though infrequent issue I have been having in QGIS. On occasion, I will begin processing a tool that will take some time to run but realize midway through that I do not want to complete the algorithm with the provided inputs. Previously I have had to end the task with my task manager since I have found that closing out of the window does not always stop the process. Ending in this way can cause data loss sometimes. Does anyone know of a way, inside the QGIS console, to stop a process?


texts:
Having issues with slope analysis. I click raster, analysis and then slope.

I enter my input and output and run the tool. The output I get doesn't appear to have been calculated correctly. Instead of seeing a continuos gradual change in slope through various colors, I just see a scattering of small squares.

It looks like only the max and min values are being displayed an nothing else (not sure).

Do I need to use the terrain analysis plugin for the tool to work? Or is there something in my workflow that I need to change/add? I am using version 3.0.3


texts: Checking if result of st_difference is empty

texts:
I am trying to select lines that completely touch a polygon (not at just a single point).
I would assume it would be a simple case of st_isempty(st_difference(geom_a,geom_b)).
However, this is not the case as you can see in the following picture.

I would like to only have the highlight yellow lines returned. What am I doing wrong?
(I could solve it with a string comparison, but I would like to avoid the extra "st_astext" invocation).
The query:
SELECT a.geometry,
       a.id as a_id,
       b.FID as b_id,
       st_astext(ST_GeometryN(st_difference(a.geometry,b.geometry),1)) as st_difference,
       st_isempty(ST_GeometryN(st_difference(a.geometry,b.geometry),1)) as st_isempty
FROM a
INNER JOIN b ON st_touches(a.geometry,b.geometry)


texts: I wonder if there is a way to interpolate xyz points in a raster within certain bounds using R? In ArcGIS this is called interpolation with barriers. I wonder if there is something similar in R? I would like to use the kriging method.


texts: Saving ST_Envelope result gives error: Geometry type (MultiPolygon) does not match column type (Polygon)

texts:
I'm trying to save the bounding box of a multipolygon into a polygon column like this:

UPDATE "mytbl" SET bb = ST_SetSRID(ST_Envelope(mp::geometry),4326) WHERE "mytbl"."id" = 1

But I receive the error:

PG::InvalidParameterValue: ERROR: Geometry type (MultiPolygon) does not match column type (Polygon)

I also tried to use the ST_ConvexHull function but it gave me the same error.


texts: Mapserver can't write in /var/www/mapserver/log/

texts:
I'm executing Mapserver 5.6.4 under Apache2 in a Ubuntu 10.04 machine.
I'm having problems writing log file. In the mapfile, I have:

CONFIG  "MS_ERRORFILE" "/var/www/mapserver/log/mylogfile.log"
DEBUG   5

Any request I make gets the same answer

msSetErrorFile(): General error message. Failed to open MS_ERRORFILE
/var/www/mapserver/log/mylogfile.log

If I change the path to /tmp/mylogfile.log it works, but then I have a
similar problem with SLD:

msHTTPExecuteRequests(): HTTP request error. HTTP: request failed with
curl error code 7 (Failed connect to mymachine:80; Operation now in
progress) for http://mymachine/sld/polygon_simplepolygon.sld
msSLDApplySLDURL: WMS server error. Could not open SLD
http://mymachine/sld/polygon_simplepolygon.sld and save it in
temporary file /var/www/mapfiles/4df0f886_22dd_0.sld.xml. Please make
sure that the sld url is valid and that imagepath and imageurl are set
properly in the map file

IMAGEPATH and IMAGEURL are properly set in Mapfile (a writable
directory and a correct url). The log directory belongs to www-data,
and the SLD directory too, and www-data has write permission in both directories.  I made a test executing, with the www-data user, a python script that creates a file in those directories and it works.

Why mapserv process can't do it? Even when I change the permissions of log directory or SLD directory to 777, mapserv can't write on them, but it can it /tmp. Why?

Many thanks in advance, and best regards,


texts: I have a set of coordinates for the four endpoints of a TIFF image. I want to add those to the image attributes so I can see the GPS location using the properties of the image.


texts: Why is data from Collections.DefaultDict(list) a string?

texts:
While creating a defaultdict(list) containing street names, x and y coordinates, my x and y entries are converted to strings. This prevents me from doing calculations on them. I've been looking around but can't find a quick fix.

Here's my code and the snippet of the data when I print it out :

import arcpy
import collections
import pandas as pd
import numpy as np
import statistics
from scipy.spatial import *

#Environments
arcpy.env.workspace = r'C:\Users\Contractuel\Documents\Etape_3'

#Local variables
fc = arcpy.env.workspace + '\SpatialJoin_Candiac.shp'
fields = ['NOMRTE','X','Y']

d = collections.defaultdict(list)

with arcpy.da.SearchCursor(fc, fields) as cursor :
    for roadname,x,y in cursor:
        d[roadname].append([x,y])

Snippet :

Rang Saint-Andre:
[['617035,993369637', '5024323,46115644'],
['616978,930218717', '5024342,95813809'],
['616935,739846513', '5024366,96582001'],
['616898,918096905', '5024394,63637209'],
['616866,027555116', '5024438,9422574'],
['616845,464712524', '5024486,52105514']]

Anybody has ever encountered this problem?


texts: Problem with Memory storage type in PyQGIS when updating with auto form

texts:
I have a layer of "Memory storage" type (Memory provider) which I
created with Python.

When editing=on and updating values using the attribute table interactively everything works as expected, i.e. the record get's updated.

Then I'm checking "Auto open form" to on. When updating the object interactively using the "Identify Features" tool the auto form opens as expected.

But after saving, I end up with two records sorted to the end, one is the unchanged original and the other record contains the changes (as well as e.g. the same id).

I didn't do anything in my code with the "Identify Features" tool. So this is a bug or at least a very strange default behaviour of PyQGIS.

Any hints?


texts: Export results of MPRI and NDVI to CSV file?

texts:
How can I export the result of MPRI or NDVI (which are rasters) to a CSV file?
I want to interpret and compare MRPI and NDVI.


texts: Showing two locations in map closer to each other than they actually are using ArcMap?

texts:
I'm working on ArcMap 10.2.2.

I need to print two maps for three municipalities, a city (Florence, Italy) and two smaller towns (Signa and Impruneta) in its vicinity. One map with the city and the other one with the two towns (it's no worth to have a map for each of these two towns).

I want the maps to have the same scale and to fit the paper and they must be comprehensive of borders polygon, basemap (only shown within the borders of the municipalities, leaving the rest of the paper blank), and point features representing my data. The way to go is to zoom to layer in map view and than set the layout the way I want, fitting the map to paper margins. I also clipped the data frame to the borders shape in order to have the basemap shown only within the borders.

The problem is that the two towns are located in different positions around the city, resulting in a large blank space left by the city in the second map (two towns map). Moreover, since I want the scale of the two maps to be the same, the two towns would partially fall outside of the paper extent. Thus, I decided to move the points and border of a town close to the other and draw a dashed line between them to make clear that's not their relative position (I also put a smaller map to show their real positions).

I know this means modifying the data (modifying the coordinates of the points etc) of the town I moved and this may not be the best way to go (I'd like advice about this), but I obtained the result I was looking for. Except for...the basemap!

Clipping the data frame doesn't not create a new clipped feature, it just hides what falls outside the clipping area. So, while I was able to move points and border, the basemap shown within the border is just the wrong area. And this is natural, since I suppose there's no way of moving pieces of basemap around, unless there's a way to export the wanted piece of basedmap as graphic or raster maybe (?!).

I hope I explained it clearily enough, however the three pictures below could help give a better idea: first picture shows map with Florence, second picture shows map with the two smaller towns (the border of Signa enclose the basemap of a wrong area) and third map shows the relative positions of the three municipality in reality.

How can I achieve what I want?


texts: How to best set up QGIS for GDA 94 zone 55 (Australia)

texts:
I am looking for help on the best way to set up QGIS for GDA 94 zone 55.
(i.e. all the relevant settings etc.)[1.8 in lieu of 2.0]

Am working off Garmin GPS/GPX data in GDA/UTM format (I access this data in QGIS usually after dumping it in Garmin Base Camp) and want to finalise maps in the GDA 94 zone 55 CRS.  I am using other ESRI/shapefile layers/data which are in GDA 94 zone 55.

Currently I am having issues with measure/buffers/grids/scales etc.
i.e. to make a 50m buffer on a layer I have to enter 50m as 0.0005 as buffer distance.
i.e. can't make grid at all in this CRS

Any help appreciated.

Please consider that the questioner [as a not for profit volunteer] does not have formal training in GIS and is currently learning QGIS "on the fly"... and so appreciates linear explanations/GUI solutions.

The questioner is extremely grateful for all the work of those contributing to QGIS/FOSSGIS/these question/answer forums etc.


texts: How to make label show up on contour vector in QGIS?

texts: I'm new to QGIS (version 2.6.0). I think my question should be easy to solve but I still get stuck at making labels show up on contour lines. These lines were generated from extracting raster (asc file). Then I tried to label those contour lines on the plot. I checked the box of "label this layer with XX" (XX is the attribute name I need) in layer properties, and clicked OK, but the there is no label showing up, I also changed the front size, but still I can't see any label. Is there anything else I need to do to make it work?


texts: Is Esri GeoPortal compatible with ArcGIS 10.2?

texts:
I am attempting to set up GeoPortal 1.2.4 and use ArcGIS 10.2 as my client desktop.  I haven't come across any documentation which says these are compatible other than the Git page.  However, I am unable to set up my software, I am getting a message which says 10.1 is not installed.

https://github.com/Esri/geoportal-server/wiki/Preinstallation-1.2.4


texts: Combine multiple Python installations coming with GIS software

texts:
On my laptop I have QGIS and ArcGIS installed. Both come with separate Python installations, and next to these two there is also a 'standard' Python installation I have for development. Now I have 3 Pythons:

Standard Python - C:\Python27
QGIS Python - C:\OSGeo4W64\apps\Python27
ArcGIS Python - C:\Python27\ArcGIS10.4

Now say I'm eg. developing a QGIS plugin and need an external library to install. It is QGIS-Python which needs that, but what if it's a common lib like matplotlib or pandas which I would like to use normally, not only in QGIS plugins? Then I must install it as well into Standard-Python. And if for ArcGIS scripting, then also for ArcGIS-Python.

I don't really want to triple my Python installations and their additional libraries. But each one mentioned here differs from others, because it has some libs necessary to run modules of said software. Therefore I can't just delete them and change Python paths in settings.

I would like to have one Python installation so I can access any module from one place. Is that possible? How to do this?

I didn't ask this question on SO because I think knowledge of GIS environment is essential here to help distinguish Pythons and know how to install them


texts: Getting centroid of closed line in MYSQL

texts:
I have a table that contains only closed lines (lines where endpoint is the same as startpoint).

I would need to find a centroid for those lines.

How would I be able to achieve it?

It does not need, to be inner centroid, it can also be just an average of all the line nodes.


texts: Efficient Outputs for GeoServer Vector Data

texts:
We have a mapping application at http://explorer.geowessex.com which uses OpenLayers and GeoServer for delivering the data. The data is potentially endless and we have very little control on what goes in. The data is served from an Oracle Database and published to GeoServer and also (in instances where it is necessary) cached off in GeoWebCache. The problem is trying to work out the correct image format to deliver the 256x256 tiles for these vector layers. We recently switched to PNG8 as this delivers tiles four times smaller than PNGs with what we thought was very little impact. Sadly we recently spotted that if we set the option Transparent to True in our GetMap request we lose any translucency that there should be on the layer, and since we are tiling layers over the top of each other, we need that option!

Here are two examples;

With 'transparent' set to false we can see the pink lines underneath the pink polygon looking thing, but the background is not transparent, so when put on a map it covers all the mapping behind it. Not good.

http://wms.geowessex.com/geoserver/wms?LAYERS=dev%3AALL_STYLE_TEST_TABLE4720,gis:ROUTE&CACHE=false&TRANSPARENT=false&FORMAT=image%2Fpng8&EXCEPTIONS=application%2Fvnd.ogc.se_xml&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A900913&BBOX=-264166.03050625,6594375.6422375,-254382.0908875,6604159.5818563&WIDTH=256&HEIGHT=256&format_options=antialias:none

With 'transparent' set to True we cannot see the pink under the pink polygon looking thing, but the background is transparent, so when put on the map you can see other layers. Better, but still not good.

http://wms.geowessex.com/geoserver/wms?LAYERS=gis:ROUTE,dev%3AALL_STYLE_TEST_TABLE4720&CACHE=false&TRANSPARENT=true&FORMAT=image%2Fpng8&EXCEPTIONS=application%2Fvnd.ogc.se_xml&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A900913&BBOX=-264166.03050625,6594375.6422375,-254382.0908875,6604159.5818563&WIDTH=256&HEIGHT=256&format_options=antialias:none

Of course, one way to solve this problem is to switch all vector data to PNGs, which work perfectly, like this

http://wms.geowessex.com/geoserver/wms?LAYERS=gis:ROUTE,dev%3AALL_STYLE_TEST_TABLE4720&CACHE=false&TRANSPARENT=true&FORMAT=image%2Fpng&EXCEPTIONS=application%2Fvnd.ogc.se_xml&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A900913&BBOX=-264166.03050625,6594375.6422375,-254382.0908875,6604159.5818563&WIDTH=256&HEIGHT=256&format_options=antialias:none

But as you will see if you look at the size of the response, this tile on it's own is 256.6KB, where as the PNG8s were 64.4KB. There is every possibility that someone might actually want 10 layers turned on at once, with each layer maybe needing about 30 tiles each depending on the users screen size you can see the amount of data that needs to be sent starts hitting many, many Megabytes a time and over the course of one session the user might have clocked a total download running into hundreds of megabytes. This is really bad and we cannot find a solution that both provides smallish tiles sizes and the ability to be transparent and translucent at the same time. Does anyone have any ideas on this? What do other mapping systems use? FYI I have looked at the all purpose GeoServer bible 'GeoServer on Steroids' as well as a number of other similar PDFs to very little avail. None seem to hit the same problem we do and I can't work out what we are doing differently! Also FYI we will be switching our OS Basemaps/Aerial Photography and other Raster datasets to JPEG soon which will save a lot of bandwidth/frustration, but we can't do the same with vectors.


texts: Populating string field based on numeric value from numeric field in ArcGIS Field Calculator?

texts:
I am new to Python and I'm trying to write a simple if then statement in Field Calculator in ArcGIS 10.2.2 that will take the values from a numeric field and add them to a string field.  I can't figure out which python function to use.  When i try this code, it works but all values in my text field become "N/A."

This is what I put in the Pre-logic Script Code:

def change(Value):
    if Value == 1:
      return "Developed"
    elif Value == 2:
      return "Agriculture"
    elif Value == 3:
      return "Forested"
    elif Value == 4:
      return "Water"
    elif Value == 5:
      return "Wetlands"
    else:
      return "N/A"


texts: QGIS working with a point and polygon layer

texts: I have 2 layers, layer A ( point data) and layer B (polygon data).   Within layer A I want a column which will indicate Yes or No if a polygon is within 500m.  I have created a 500m buffer for layer b and have currently used the geomtouchs function to select those polygons.  I can't get it to fill the column with Yes or No though.  Has anyone any other methods they would use?


texts: How to keep native spatial resolution while making a chart?

texts:
I have difficulties on producing chart in native spatial resolution. Please, look at the below script. The script allows to transform hourly GSMaP data into daily data and extract the values for given period as chart. But, to produce the chart I am required to specify the spatial scale which is "10000" in my script, while I would like to extract the data as chart in native resolution. If I delete the "10000" from chart line and run the script, GEE shows error:

Please, help me to improve the script to be able to produce chart in native resolution.

var pts = /* color: #d63000 */ee.Geometry.Point([65.945, 42.027]),
    aoi = ee.Geometry.Polygon(
        [[[64.62748181270217, 43.19201550614239],
          [64.62748181270217, 40.77519027755922],
          [67.43998181270217, 40.77519027755922],
          [67.43998181270217, 43.19201550614239]]], null, false);

// Start and End Dates
var inidate = ee.Date('2014-02-01')
var enddate = ee.Date('2014-03-30')

//Import GSMaP data
var img_1 = ee.ImageCollection("JAXA/GPM_L3/GSMaP/v6/reanalysis")
    .filterDate('2000-03-01', '2014-03-01').select('hourlyPrecipRateGC');

var img_2 = ee.ImageCollection("JAXA/GPM_L3/GSMaP/v6/operational")
    .filterDate('2014-03-01', '2030-12-31').select('hourlyPrecipRateGC');

var imgcol = img_1.merge(img_2)
    .filterDate(inidate, enddate);

// Difference between start and end in days
var difdate = enddate.difference(inidate, 'day')

// Time lapse
var lapse = ee.List.sequence(0, difdate.subtract(1))
var listdates = lapse.map(function(day){
  return inidate.advance(day, 'day')
})

// Function for producing daily image collection
var daily_gsmap = ee.ImageCollection.fromImages(listdates.map(function(day){
  day = ee.Date(day);
  var daily_collection = imgcol.select('hourlyPrecipRateGC').filter(ee.Filter.date(day, day.advance(1, 'day')));
  return daily_collection.sum().clip(aoi).rename('daily_rainfall').set('date', day.format('YYYY-MM-dd'))
        .set('system:time_start', day.millis()).set('system:time_end', day.millis());
}));

//Making chart and printing it
var chart_daily = ui.Chart.image.seriesByRegion(daily_gsmap, pts, ee.Reducer.first(), 'daily_rainfall', 10000)
  .setChartType('ColumnChart');
  chart_daily.setOptions({
    title: 'Daily rainfall',
    vAxis: {title: 'Rainfall (mm)'},
    hAxis: {title: 'Date', format: 'yyyy-MM-dd'},
  });

print(chart_daily);


texts: How to open a file geodatabase with tables and relationships in QGIS?

texts: How can I open or connect to an ESRI file geodatabase in QGIS? I don't want to only see the features (shapes), I want to make sure that tables and relationships are opened correctly, too.


texts: Rasterio and WarpedVRT to go: Geo > UTM (compute stuff) > GEO without changing final extent?

texts:
I try to take 1arc (1 arc-second res.) SRTM1 elevation data, project it into UTM, do slope, aspect and other calculations that require a meter-based projection, and finally reproject the results back to the original projection. At the moment I use a chain of gdal commands at the bash and a set of python files to achieve this, but is's messy and takes a lot of disk space.

So I plan to do all steps in a single script using rasterio. My basic idea was:
- open LatLon source file with WarpedVRT, converting it to UTM on the fly
- create multiple new raster bands (slope, aspect, other) that I add to the WarpedVRT instance
- open a second WarpedVRT instance (this time projecting the extended WarpedVRT_UTM) to get my new data back to LatLon and the original resolution and extent

But by checking the array shapes I see the following:
- original GeoTiff (LatLon): 3601x3601 px
- WarpedVRT (UTM): 3904x3435 px [this is OK an expected]
- WarpedVRT (LatLon): 3746x3788 px [this should be the original res and location]

I tried to get some inspiration from the commandline tools (basically window clip the 2nd WarpedVRT with something that resembles the '--like file' functionality of rio.clip but I'm stuck.

Any points if I am thinking along the right lines here? Or should I use another class?

My first attempts were along the following lines:

dst_crs = "EPSG:%d" % calc_utm_epsg(coord)

with rasterio.open('s30_w072_1arc_v3.tif') as src:
    src_crs = src.get_crs()

    print src.shape
    # shape (3601, 3601)

    with WarpedVRT(src, dst_crs=dst_crs, resampling=Resampling.bilinear) as vrt_utm:
        print vrt_utm.shape
        # shape (3904, 3435)

        # (1) clone vrt for slope, aspect to hold results
        #     (copy profile etc.),
        #     maybe place new data bands into one vrt
        #
        # (2) do slope, aspect calulations
        # (...)

        # back-warp cloned vrt to original shape
        # change vrt_utm to slope_aspect_vrt

        with WarpedVRT(vrt_utm, dst_crs=src_crs, resampling=Resampling.bilinear) as vrt_src:
            print vrt_src.shape
            # shape (3746, 3788)

I'm using RasterIO v 1.0a9, GDAL2 on a Mac (homebrew)


texts: Starting ArcMap says system clock has been set back more than 24 hours?

texts:
In order to test my log4net logging daily, I changed my clock to be 2 days later.

When I changed back the clock to normal, I cannot start ArcMap, I get a message:

  The system clock has been set back more than 24 hours

I found a possible solution provided by Esri, but it not work.

Is there another solution for it?

I am using Wndows 7 and ArcGIS 10.2 for Desktop.


texts:
I am using Osm2Po to find routes between two points.

I am trying to find routes based on the difficulty specified by the user.

Basically, I would like to calculate the cost based on some parameter in the way. For now this parameter could be the type of road e.g. if its a primary road the cost would increase. If it was a secondary road it the cost would decrease especially. Any road with a cyclelane would have the best cost.

I have been looking into overriding the DefaultRouter methods: traverse, findpath and calcEdgeCost.

I was thinking that I could get the clazz type from the edge being currently examined and set its cost using calcEdgeCost.

One problem is that for some reason calling getVertexClazzes returns nothing but zeroes. getEdgeFlags returns null.

Here is part of my osm2po.config which might be relevant

#wayTagResolver.class = com.to.another.package.YourWayTagResolver
    #nodeTagResolver.class = com.to.another.package.YourNodeTagResolver
    #relationTagResolver.class = com.to.another.package.YourRelationTagResolver

# Properties of the DefaultWayTagResolver use this prefix.
wayTagResolver.prefix = wtr

# Custom flags with ascending binary values 1, 2, 4, 8 ...
# You may define up to 32 Flags (Bits).
wtr.flagList = car, bike, foot, rail, ferry

# final decision; only allow ways with these flags
wtr.finalMask = bike

# very special hint for level_crossing modification
wtr.shuttleTrainMask = rail|car

# Main-Tag definitions. Params 1-4:
# 1) concurrent order
# 2) class (1-127)
# 3) default speed in kmh
# 4) allowed transportation type (optional) - since v4.5.30

wtr.tag.highway.motorway =       1, 11, 120, car
wtr.tag.highway.motorway_link =  1, 12, 30,  car
wtr.tag.highway.trunk =          1, 13, 90,  car
wtr.tag.highway.trunk_link =     1, 14, 30,  car
wtr.tag.highway.primary =        1, 15, 15,  car|bike
wtr.tag.highway.primary_link =   1, 16, 30,  car
wtr.tag.highway.secondary =      1, 21, 15,  car|bike
wtr.tag.highway.secondary_link = 1, 22, 15,  car|bike
wtr.tag.highway.tertiary =       1, 31, 12,  car|bike
wtr.tag.highway.residential =    1, 32, 12,  car|bike
wtr.tag.highway.road =           1, 41, 14,  car|bike
wtr.tag.highway.unclassified =   1, 42, 12,  car|bike
wtr.tag.highway.service =       1, 51, 5,   car|bike
wtr.tag.highway.living_street = 1, 63, 7,   car|bike|foot
wtr.tag.highway.pedestrian =    1, 62, 5,   bike|foot
wtr.tag.highway.track =         1, 71, 10,  bike|foot
wtr.tag.highway.path =          1, 72, 10,  bike|foot
wtr.tag.highway.cycleway =      1, 81, 15,  bike
#wtr.tag.highway.footway =       1, 91, 5,   foot
#wtr.tag.highway.steps =         1, 92, 5,   foot
wtr.tag.route.ferry =            2, 01, 10,  ferry
wtr.tag.route.shuttle_train =    2, 02, 50,  rail|car
wtr.tag.railway.rail =           3, 03, 50,  rail


texts: Using MOD16A2 in GEE to extract best quality of PET

texts:
I am beginner with Google Earth Engine programming. I want to know how I can use MOD16A2 to extract good quality values of potential evapotranspiration data (MOD16A2)? because I extract the PET of MOD16A2 and I compared with PET calculated by FAO56 PM, I find the bug difference between this parameters. My study area is agricultural area.
This is my code :
//Load the dataset
var dataset = ee.ImageCollection("MODIS/006/MOD16A2")
               .filter(ee.Filter.calendarRange(2001,2020, 'year'))
               .filterBounds(aoi);// aoi= study area

var PET = dataset.map(function(img){
  var PETmm = img.select('PET')
                 .multiply(0.1);
  return PETmm;
});
print('8day potentiel evapotranspiration in image collection',PET);
print(PET.size());
print(aoi.size());
print(PET.toBands());
// reduce the ImageCollection to singleimage with  multiple band
function imageCollectionToBands(imageCollection) {

  // put all images in 1 image with multiple bband
  var img = imageCollection.toBands();

  // Rename bandnames
  var names = img.bandNames();
  var newNames = names.map(function(name){
    var ind = names.indexOf(name);
    return ee.String(names.get(ind)).slice(0,10);
  });
  return img.rename(newNames);
}

var imagePet= imageCollectionToBands(PET)
                   .where(imageCollectionToBands(PET).gt(100),ee.Image.constant(0));
print('8day PET', imagePet);
var marocPet = imagePet
  .reduceRegions({
    collection: aoi.select(['Nom_Locali']),
    reducer: ee.Reducer.mean(),
    scale: 1000
  });

print('Feature collection of PET', marocPet);
// Save Value in google_drive
Export.table.toDrive({
  collection: marocPet.select(['.*'],null,false),
  description: 'PETStations',
  fileNamePrefix: 'PET_Stations',
  folder : "GoogleEarthEngine Data",
  fileFormat: 'CSV'
});


texts: Aligning points in WGS84 to centroid of features in projected CRS (EPSG:28992) using QGIS

texts:
I have a set of points in WGS84 representing the rainfall over the Netherlands, as acquired by radar towers. These points represent the centroid of a pixel in WGS84 and they are separately by (roughly) 1km.

In the project I am working, I would like to use a projected CRS, in this case EPSG:28992 (official Dutch CRS). I created a 1km grid in EPSG:28992 and I overlaid the points, after transforming them from WGS84 to EPSG:28992. You can see an image illustrating the problem below.

As expected, the centroid of the reprojected WGS84 points does not coincide with the centroid of the pixels in EPSG:28992. Hence, I believe the task to do here is to align each of these points with the closest centroid of the pixels in EPSG:28992. However, I am not sure whether this reasoning is correct, because since there is this inclination of the points I might assign a red dot to the incorrect pixel.

Can you guide me on how to do this properly using QGIS?

EDIT: The purpose of aligning is because the rest of the analysis will be done using numpy arrays (matrices), and I think it is easier to prepare the data in a projected CRS, since then each cell of the array corresponds to each pixel of the CRS. The overlaid grid is just to illustrate what the problem is.


texts: Calculating boundary around all land on Earth?

texts:
I hope this is right here in Earth Science, it's a bit of a Computer Science-y problem, too:

I want draw (or really just calculate the total area) of a boundary of X miles (with X varying between around 200-400) around all of Earth's landmass. It would be very much like the problem of drawing Territorial waters, only that I don't care about countries and nations at all. I just want a line around all landmass, including islands etc and then calculate the area of that.

Anyone knows a good way how to go about this?


texts: Can't install Google Earth Engine plugin QGIS 3.16?

texts:
When I try to install the Google Earth Engine plugin in QGIS 3.16 I have this message:
"Couldn't load plugin 'ee_plugin' due to an error when calling its classFactory() method "


texts: Trigger a layers 'on hover' effect through another layer in Leaflet

texts:
Im trying to make a map containing two separated basemaps one with terrain and second with labels, and with vector data layers in between. The problem is when I try use hover/click events on data layers because the label layer on top is blocking them.
Is there a way to ignore specific layer in mouse events?

the code based mostly on choropleth map leaflet tutorial:

var baseMap1 = L.esri.basemapLayer("Oceans");               //terrain without labels

var baseMap2 = L.tileLayer('https://api.mapbox.com/..., {   //mapbox labels only
               pane: 'shadowPane'});                        //just for testing , it will have its own pane on top of everything else

var map = L.map('map',    {
    layers: [baseMap1, baseMap2],
    });

var vectorLayer = L.geoJson(vectors,  {style: style, onEachFeature: onEachFeature,}).addTo(map);

function onEachFeature(feature, layer) {
    layer.on({
    mouseover: highlightFeature,
    mouseout: resetHighlight,
    });
}

function highlightFeature(e) {
        this.setStyle({
        weight: 5,
        color: '#666',
        dashArray: '',
        fillOpacity: 0.7
        });
       }

 function resetHighlight(e) {
        vectorLayer.resetStyle(this);
        this.bringToBack();
        }

I figured there is an 'interactive' option for vector layers with polygons that can do what I need (polygons are "transparent" for mouse events) but it doesn't work for tileLayer

JS fiddle - http://jsfiddle.net/qkyejuxy/
comment 4th line and set the labels below vector layer  and the hover works but labels are not visible


texts: How to check if a point is near or inside polygon?

texts:
I am looking for a simple library which can check, whether a point is inside a polygon or not. Earlier I found Java Spatial Index library, but couldn't figure out how it works.

I also tried openmap, again without success.

I am totally new to GIS. My application has polygons. I need to check if my point is in a polygon or new polygon(s).

In open map I wrote as follows:

package abc.poi;

import com.bbn.openmap.geo.Geo;
import com.bbn.openmap.geo.OMGeo;

public class PoiTest {

    public static void main(String[] arg) {
        Geo geo1 = new Geo(1,1);
        Geo geo2 = new Geo(3,2);
        Geo geo3 = new Geo(4,4);
        Geo geo4 = new Geo(3,4);
        Geo geo5 = new Geo(4,3);
        Geo geo6 = new Geo(1,1);

        Geo listGeo[] = {geo1,geo2,geo3,geo4,geo5,geo6};

        OMGeo.Polygon poly = new OMGeo.Polygon(listGeo);

        System.out.println(poly.isPointInside(new Geo(2,2)));
    }
}

It always returns false. What am I missing?


texts:
I am trying to save this print composer as PNG

When I open the saved file, this is the result

The same thing happens if I save to pdf. Any idea what's going on?


texts: Does this routable network include oneways information?

texts:
I successfully imported osm road data using osm2po. I visualized the routable road network using Qgis. The table of the layer includes the following attributes: id|osm|name|osm_source|osm_target|clazz|flags|source|target|km|kmh|cost|reverse_cost|x1|y1|x2|y2
Does this routable network include oneways information or do I need to further process the data in order to get this information?


texts:
I am trying to export a csv of daily mean cloud cover in the Netherlands according to S5P data. I have calculated the mean cloud fraction, but the exported csv shows many, many rows - none of them are cloud_fraction (image below). What am I missing?
//  1. Setup
var countries = ee.FeatureCollection("FAO/GAUL/2015/level0");
var NL = countries.filterMetadata('ADM0_NAME', 'equals', 'Netherlands');
var start_period = ee.Date('2018-07-01');
var end_period = ee.Date(new Date().getTime());

var clouds = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_CLOUD')
  .select('cloud_fraction')
  .filterDate(start_period, end_period);

var cloudNL=clouds.map(function(im){
   return im.clip(NL)//.resample('bilinear').reproject('EPSG:4326', null, 2000);
});

//  2. Mean cloud cover dictionary for entire collection
var cloudlist = cloudNL.map(function(image){
  var meandict = image.select('cloud_fraction').reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: NL
  })
  return image.set(meandict);
});

Export.table.toDrive(cloudlist.sort('cloud_fraction'), 'NL_mean_cloud_fraction');


texts: How to install pgRouting in Ubuntu 14.04?

texts: I want to install pgRouting in Ubuntu 14.04 LTS. I added the respective repository (https://launchpad.net/~georepublic/+archive/pgrouting-unstable?field.series_filter=trusty) but failed to install the package in the terminal so far. Did someone encounter this before? Thanks for every piece of advice.


texts: Converting map extent from 3857 projection to 26331 projection

texts: I am trying to convert extent which is in 3857 projection into 26331 projection in OpenLayers by ol.proj.transformExtent(map.getView().calculateExtent(map.getSize()),'EPSG:3857','EPSG:26331')


texts: No data showing when adding delimited text layer with lat/long

texts:
See screenshot of Data Source Manager of Delimited Text input. Sample data look fine.

When I press 'Add' no data shows. World Map layer is EPSG-3857-WGS 84 / Pseudo Mercator. Am I using the wrong CRS? Using MacOS 10.14.6 Mojave and QGIS 3.12.1.


texts: Tool in ArcMap for capturing ROIs of same size?

texts:
I'm trying to run supervised classifications and need all ROIs to be the same size. I'd like to capture the ROIs in ArcMap and so far can only find a tool which allows you to draw your ROIs by clicking and dragging to create polygons.

In ENVI and Erdas Imagine, there is a tool for capturing ROIs to ensure they are all the same size. For example if you want each ROI to be 9 pixels, when you're creating signature files, the cursor is 3x3 pixels. Then where you click on your raster, the 9 pixels are automatically captured without the need to draw a polygon. I can't remember what it's called as it's been a while since I used it.

After capturing the ROIs, I plan to export them as a tiff to use in EnMap

Does anyone know if such a tool exists in ArcMap?


texts: r.resamp.interp is not giving correct results

texts:
I am doing resampling uisng r.resamp.interp in GRASS GIS, but, after applying each time It reduces the size of raster. I am using this command
g.region raster=raster
r.resamp.interp input=raster output=resampled_raster method=bilinear

Original raster

Resampled rater


texts: How to identify landuse change with Sentinel data - specifically grassland/woodland to arable cropland changes

texts:
I have an assignment to assess whether it is viable to use sentinel data to monitor  illegal changes of landuse from Grassland/Woodland to arable cropland.

This is the first time my environmental agency has shown interest in Sentinel open-data and I would like to show that it is cost effective in financial terms as well as in terms of human resources. The group consists of a GIS manager (myself) and 2 ecologists. I would also like to show whether the results are valuable in practical terms (i.e reliable, so that we can use the results to compliment ground observations to indicate whether landowners are illegally turning grassland into arable land, or too complex to carry out without external assistance.

I have been to a 3 day Copernicus seminar and know the basics of snap, and also have good ArcGIS and FME experience (SAGA, and QGIS as well). I have only minimal knowledge of red-edge, infrared analysis.

Does anyone have any steps or information for beginners? I can only find pretty complicated material which is rather too detailed.


texts: Color Ramp of Raster Rendering in MapServer

texts:
Does this proposed new syntax work yet? http://mapserver.org/development/rfc/ms-rfc-6.html#proposed-new-syntax

Currently, I am still using this method for rendering the raster in Mapfile.

STYLE
   COLORRANGE 0 0 0  255 255 0 # black to yellow
   DATARANGE 0.0 100.0
   RANGEITEM "foobar"
END

However it doesn't work so well since it is linear interpolation and equal interval. In reality, the histogram of raster shows that data doesn't distribute with equal interval. In ArcGIS, the default method for classified raster rendering is Natural Breaks (Jenks).

Is there any new improvement for that in MapServer? Thanks


texts: Plotting irregular raster object with tmap

texts:
I would like to use tmap R package to plot a RasterStack object. The problem is that the coordinates are not real coordinates, but rather climate variables, so I can not specify a valid reference coordinate system. I can plot this object in ggplot and lattice, but I would really like to do it in tmap.
> temp_object
class      : RasterStack
dimensions : 400, 550, 220000, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.1, 10  (x, y)
extent     : -20, 35, 0, 4000  (xmin, xmax, ymin, ymax)
crs        : NA
names      :          X1
min values : 0.001695915
max values :   0.1038332

> tm_shape(temp_object, bbox=tmaptools::bb(matrix(c(-20,0,35,4000),2,2))) + tm_raster()
Error: The projection of the shape object temp_shp is not known, while it seems to be projected.

If I try to specify a projection, I do get my plot, but the aspect of x and y is not OK and I did not find any argument in tmap R package, which could give me a desired aspect ratio.
crs(temp_object) <-  "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
tm_shape(temp_object, bbox=tmaptools::bb(matrix(c(-20,0,35,4000),2,2))) +
  tm_raster(legend.show = FALSE) + tm_grid()

Is there any argument in tmap, which would enable me to plot this raster with random aspect ratio? An example of such argument in ggR R package is argument coord_equal = FALSE.
  library(RStoolbox)
  ggR(temp_object,  coord_equal = F, geom_raster = 1)


texts:
I'm using Mapnik to show a map with hillshading. It works well when displaying it up to my tiff's native resolution and a little bit above, for instance at zoom 13 (dataset is SRTM30):

But then if I continue to zoom in, the hillshade isn't smooth anymore and looks like vector polygons, like on this screenshot at zoom 16:

I'm currently using bicubic scaling but I tried every single other option without seeing much improvement. Upscaling my raw image at the same resolution with Photoshop works perfectly, with smooth fading.

Is there a way to improve this result? I tried to add agg-stack-blur without success.