texts: Render WMS at fixed scale (static image effect)

texts: I have a point layer that I'm rendering as a WMS using GeoServer, the points look like a continuous surface and 'blend' great at about 1:10 000. I want to keep the WMS rendered at that scale, irrespective of the client's zoom level. So if the user is zoomed out at 1:5 000 000 or zoomed in at 1:500, the WMS must remain at zoom level 1:10 000. I currently deal with this using scale ranges and rules in an SLD, but surely there must be a better way?


texts: Allow CORS on GeoServer with Tomcat and Nginx

texts:
I've used this procedure for install GeoServer on a server based on Ubuntu 20.04 with Tomcat and Nginx. At the end I've putted under HTTPS my subdomain.
I can see without problems the WMS services and infact I've created this map that use eight WMS link.
The problem comes when I try to pick up the informations from the pixels because I see into console log many problems about CORS origin.
So, I've followed this procedure to activate CORS on GeoServer. After the edit of web.xml I've reloaded GeoServer clicked on Reload into web user interface but GeoServer goes offline with the error:

HTTP Status 404 – Not Found
Type Status Report
Description The origin server did not find a current representation
for the target resource or is not willing to disclose that one exists.
Apache Tomcat/9.0.37

I've restarted Tomcat and Nginx with:
systemctl restart tomcat
systemctl restart nginx

and now I see another error:

HTTP Status 404 – Not Found
Type Status Report
Message The requested resource [/geoserver/] is not available
Description The origin server did not find a current representation
for the target resource or is not willing to disclose that one exists.
Apache Tomcat/9.0.37

How I can solve this problem?
EDIT 1:
Inside the logs (geoserver/data/logs) there aren't indication about the problem. If I comment again the rows edited inside web.xml and use systemctl restart tomcat it is possible to use again GeoServer.
EIDT 2:
Inside web.xml I've this rows that are commented:
   <!-- Uncomment following filter to enable CORS
<filter>
        <filter-name>cross-origin</filter-name>
        <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>
       <init-param>
           <param-name>chainPreflight</param-name>
           <param-value>false</param-value>
       </init-param>
       <init-param>
           <param-name>allowedOrigins</param-name>
           <param-value>*</param-value>
       </init-param>
       <init-param>
           <param-name>allowedMethods</param-name>
           <param-value>GET,POST,PUT,DELETE,HEAD,OPTIONS</param-value>
       </init-param>
       <init-param>
           <param-name>allowedHeaders</param-name>
           <param-value>*</param-value>
       </init-param>
    </filter>
-->
.
.
.
   <!-- Uncomment following filter to enable CORS
    <filter-mapping>
        <filter-name>cross-origin</filter-name>
        <url-pattern>/*</url-pattern>
    </filter-mapping>
-->

I uncomment the rows and then save the file.


texts: how to Zoom to max extent in openlayers3?

texts:
How to zoom to max extent on click of button in openlayers-3?

I tried with

function zoomtoextent(){
//map.zoomToMaxExtent();
if(globalbBox==null)
    {
    var bBox = new OpenLayers.Bounds(
           68.14806781208101, 20.118149140613397,
        74.47891169978561, 24.722655175866358
    );
    map.zoomToExtent(bBox);
    map.zoomTo("7");
    globalbBox = bBox;

    }
else
    {
        map.zoomToExtent(globalbBox);
    }
}

This was working fine with openlayers-2 but with openlayers-3 it is not working.


texts:
I'm working with a dataset that contains millions of polygons representing leased properties. In many cases, there are multiple polygons that are drawn on top of one another (representing multiple leases for the same property).

Currently I'm working on a process that will group these so that only one polygon geometry is used to represent multiple lease counterparts (same property leased multiple times). I have come up with a python script that uses an arcpy Select by Location, and compares the XY centroid coordinates, perimeter, area, and extent of overlapping polygons to apply a unique ID to these groups (basically). This works just fine, but it takes forever.

I'm wondering: would it be possible to use SQL (Using SQL Server Management Studio) to identify overlapping polygons (based on buffer of centroid) and group by other spatial attributes of the input polygon?

I have made an attempt at getting this to work, but got too many groups. I have around 41,000 polygons (for the first county) and should have fewer than 14,000 'unique' polygons or groups.
I have access to FME and ArcGIS Advanced Licenses, so can anyone provide an alternative to using an arcpy.da Search cursor to select by location and compare selected polygon attributes in python?

Here is my first attempt at the SQL Query. The projected coordinate system uses meters, thus the +- 11 in the query for the extents.

    SELECT a.Poly_ID
    FROM
    (
    SELECT [OBJECTID]
          ,[SHAPE]
  ,[calc_acres]
  ,[perimeter]
  ,[CentX]
  ,[CentY]
  ,[X_Max]
  ,[Y_Max]
  ,[X_Min]
  ,[Y_Min]
  ,[Poly_ID]
  ,[leaseId]
  ,geometry::STGeomFromText('POINT (' + convert(varchar,[CentX]) + ' '+ convert(varchar,[CentY]) + ')', 5070) as cent_geom
   FROM [poly_topology].[gis_mfg].[LTPOLYGONS]) a
  left join (
SELECT [OBJECTID]
  ,[SHAPE]
  ,[calc_acres]
  ,[perimeter]
  ,[CentX]
  ,[CentY]
  ,[X_Max]
  ,[Y_Max]
  ,[X_Min]
  ,[Y_Min]
  ,[Poly_ID]
  ,[leaseId]
  ,geometry::STGeomFromText('POINT (' + convert(varchar,[CentX]) + ' '+ convert(varchar,[CentY]) + ')', 5070) as cent_geom
 FROM [poly_topology].[gis_mfg].[LTPOLYGONS]) b on
  a.cent_geom.STBuffer(6).STContains(b.cent_geom) = 1 and
  b.perimeter > (a.perimeter - (a.perimeter * 0.03)) and b.perimeter < (a.perimeter + (a.perimeter * 0.03)) and
  b.calc_acres > (a.calc_acres - (a.calc_acres * 0.03)) and b.calc_acres < (a.calc_acres + ( a.calc_acres * 0.03)) and
  b.X_Max > (a.X_Max - 11) and b.X_Max < (a.X_Max + 11) and
  b.Y_Max > (a.Y_Max - 11) and b.Y_Max < (a.Y_Max + 11) and
  b.X_Min > (a.X_Min - 11) and b.X_Min < (a.X_Min + 11) and
  b.Y_Min > (a.Y_Min - 11) and b.Y_Min < (a.Y_Min + 11) and b.leaseId != a.leaseId
  group by a.Poly_ID
  order by a.Poly_ID`

UPDATE:
After implementing FelixIP's idea into my process, I found the following link: http://clear.uconn.edu/tools/Shape_Metrics/method.htm

This Shape Metrics tool is old, but the methodologies seem more robust than relying on perimeter, area, extent, and centroid alone. These values are susceptible to outliers, so not every polygon in a group is caught (mainly resulting from extraneous boundaries where two polygons were dissolved but the boundaries were not completely overlapping).


texts:
I'm looking to find commercially available admin boundaries (0 through lowest available) for the world that I can place into a user-friendly interface that allows the non-technical user (ie someone who has never heard of a shapefile) to select what boundaries are preferred for generating the map visuals.

Could someone please point me in the direction of what vendors offer this kind of dataset that could be pulled into my own commercial application? Is this a Tomnod kind of thing? Google? ESRI? My assumption is that OSM availability is too patchy for this.


texts:
I have a line made of 17959 nodes. I want to delete several range of nodes from that line with a python command. For example I want to delete all nodes between node 5530 to node 5861 and connect node 5530 to node 5861. Then I want to delete all nodes between node 14094 to node 14398 and connect node 14094 to node 14398. Here is the raw code I am using now.

layer = iface.addVectorLayer("input.shp","","ogr")
layer.startEditing()
QgsVectorLayerEditUtils::deleteVertex(QgsFeatureId "input.shp", int vertex)

I am using QGIS 2.18.15


texts: Python scripts execution in ArcGIS Runtime WPF app

texts:
DEAR PEOPLE FROM THE FUTURE: Here's what we've figured out so far…

I am facing a quite interesting requirement:

  "run python scripts within an ArcGIS Runtime app".

checking the concepts in the ArcGIS Resources for ArcGIS Runtime SDK for WPF i can read:

  Geoprocessing in the runtime provides you with a framework for
  answering your GIS questions. The API works with geoprocessing though
  geoprocessing service tasks. These tasks are either hosted online,
  through an ArcGIS Server, or can be hosted locally, by providing the
  local server with a runtime enabled geoprocessing package (GPK).

  Geoprocessing services can host one or more tasks, which come from a
  tool, script, or model, authored in ArcGIS for Desktop and then
  published to ArcGIS for Server or packaged as a geoprocessing package

but also I can read the following statement

  Not all geoprocessing tools or functionality are supported in the
  runtime local server. Please see the supported geoprocessing tools
  list.

a few lines after is the next statement:

  Python scripts   Python is an open-source, general-purpose programming
  language used as a scripting language in ArcGIS geoprocessing. Python
  scripts can be shared for use with the runtime as geoprocessing
  packages. Scripts can be turned into script tools, and then packaged
  into a geoprocessing package. The ArcGIS Runtime contains its own
  imbedded Python

I understand that I can run python scripts using the geoprocessing packages, but not all the functionality that I can get from ArcPy library is available, Am I wrong?

is there a known workaround that allows to run python scripts? (with out using ArcObjects, just ArcGIS Runtime)


texts:
I have updated the x,y,z coordinates for my LAS file, using the solution from here, however the modified date for the file does not change, it stays the same.
Is it possible to update it?


texts:
I am trying to insert spatial data from GeoJSON to PostgreSQL database and before it I am to reproject geom from EPSG:4326 to EPSG:102027. When working on Windows there were no problems with CRS, but now I get in Python

  ERROR 6: EPSG PCS/GCS code 102027 not found in EPSG support files.  Is
  this a valid EPSG coordinate system? ERROR 6: No translation for an
  empty SRS to PROJ.4 format is known.

Probably the error connected with:

That I moved to Ubuntu and maybe I missed some packages?
The venv (on Windows I didn't use it)

I have installed pygdal 2.2.2.3 package, GDAL 2.2.2

in database I have POSTGIS=2.1.8

  EPSG:102027 is an 'Asia_North_Lambert_Conformal_Conic' projection with
  description (I have it in QGIS!) +proj=lcc +lat_1=15 +lat_2=65
  +lat_0=30 +lon_0=95 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

I found the same problem but no solution, this issue offers to update EPSG database - I already have GDAL 2.2.2...

Strangely, I really haven't srid = 102027 in spatial_ref_sys table!

I am puzzled, maybe any ideas?


texts: GDAL: wrong corner coordinates from KNMI grib file

texts:
I'm trying to parse KNMI grib files to extract weather forecasts for The Netherlands. I use gdal.netcore which is a .NET wrapper for gdal 3.2.0. On the website of KNMI they say the following:

Rekenrooster: Lambertprojectie
Product resolutie: 0,037° west-oost en 0,023° noord-zuid (300x300 punten)
Grid representatie: Regulier Lat/Lon
Gebiedsdefinitie: noord 55,877° N.B. zuid 49,0° N.B. / oost 11,063° O.L. en west 0,0° W.L.

However when I run
var info = Gdal.GDALInfo(dataset, new GDALInfoOptions(null));
Console.WriteLine(info);

I get the following information:
Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: /home/HA40_N25_202103190600_00000_GB
Size is 300, 300
Coordinate System is:
GEOGCRS["Coordinate System imported from GRIB file",
    DATUM["unnamed",
        ELLIPSOID["Sphere",6367470,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (-0.018500000000000,55.888500000000001)
Pixel Size = (0.037000000000000,-0.023000000000000)
Corner Coordinates:
Upper Left  (  -0.0185000,  55.8885000) (  0d 1' 6.60"W, 55d53'18.60"N)
Lower Left  (  -0.0185000,  48.9885000) (  0d 1' 6.60"W, 48d59'18.60"N)
Upper Right (  11.0815000,  55.8885000) ( 11d 4'53.40"E, 55d53'18.60"N)
Lower Right (  11.0815000,  48.9885000) ( 11d 4'53.40"E, 48d59'18.60"N)
Center      (   5.5315000,  52.4385000) (  5d31'53.40"E, 52d26'18.60"N)

This is quite close (55.8885 vs 55.877) but nonetheless inaccurate. Even if I try to project it, it still shows the same coordinates.
GdalBase.ConfigureAll();
using var dataset = Gdal.Open(@"D:\HA40_N25_202103190600_00000_GB", Access.GA_ReadOnly);
var info = Gdal.GDALInfo(dataset, new GDALInfoOptions(null));

Console.WriteLine(info);

double[] adfGeoTransform = new double[6];
dataset.GetGeoTransform(adfGeoTransform);

var upperX = adfGeoTransform[0];
var upperY = adfGeoTransform[3];

var lowerX = upperX + (dataset.RasterXSize * adfGeoTransform[1]);
var lowerY = upperY + (dataset.RasterYSize * adfGeoTransform[5]);

var projection = dataset.GetProjection();
var source = new OSGeo.OSR.SpatialReference(null);
source.ImportFromWkt(ref projection);

var target = new OSGeo.OSR.SpatialReference(null);
target.ImportFromEPSG(4326);

var transform = new OSGeo.OSR.CoordinateTransformation(source, target);
double[] data = new double[3];
transform.TransformPoint(data, upperX, upperY, 0);

Console.WriteLine(data[0].ToString() + " " + data[1].ToString());

So what am I doing wrong here?


texts: Using GRASS algorithms in Sextante under ArcGIS Desktop?

texts:
Sextante for ArcGIS comes with SAGA modules.

I know that I can use GRASS algorithms in Sextante under gvSIG.

But how can I use GRASS algorithms in Sextante under ArcGIS Desktop?

I tried to look for any instructions, but could not find any. Any suggestions?


texts: How to set a default value for attributes in shapefiles in QGIS?

texts:
In QGIS 2.14.6 I would like to define a default value for an attribute in a shapefile (or spatialite layer as well) which is applied if the user doesn't explicitely sets the value.

Is there a way to achieve this using standard tools/setting (without python)?


texts: Dynamic text in print composer from layer attributes

texts:
I have prepared my page layout in print composer, and in my title i want to generate a dynamic texts which is related to an attribute layer which i have filtered previously.

For example:

I have layer1(this layer has only one row after filtering) with attributes att1,att2,att3
I want a dynamic in my composer that puts a text label att1-att2 or just only att1

This is a screenshot of my layer and their attributes which from i want to make dynamic text.


texts:
I can not seem to create a relationship class in ArcSDE, or at least one that works.

I am succseful in creating the class however ArcMap does not recongnize it (there are no related tables shown) and if I click on relationship classes properties in the context menu in ArcCatalog I get the following error message:

Again the relationship class is created succesfully the results are included. I get the same if I use the catalog tree, toolbox or python script.


texts: Creating DSM and DTM from a LAS using a freeware solution

texts:
Is it possible to export a surface and terrain model files from an aerial, LiDAR point cloud (LAS) set?

I am using SAGA GIS for very basic operations, but I couldn't find this option (and I have a hunch SAGA GIS should support both operations).


texts:
I would like to know, whether and how it is possible to import only certain layers of a .dxf file. I made a code which imports all points and lines. However, this code can only import everything at once.
At the moment, I simply delete everything unwanted from my .dwg AutoCAD file before I export it. However if its possible to only import certain layers like for exaple 'Boreholes' and 'Anomalies', it would greatly increase efficiency.
I use QGIS version 2.8.8 'Wien'


texts: Error code 429 and Google Maps Satellite Imagery blank screen on SASPlanet

texts:
I am using SASPlanet (Version: Nightly 190426) to download satellite imagery from google maps so I can save them locally and use them offline. However, after several tiles of downloading and stitching, the images refuses to show, leaving only blank areas with "Download error: unknown status code 429" floating around.

This was also affects maps.google accessed from chrome: only roads and annotations appears when switching from map view to satellite view but no imagery. Clearing chrome cookies solved the problem on chrome, but not on SAS.

What seems to be the cause and how to fix it?


texts: Using QGIS polygonize tool to create discreet slope polygons

texts:
I am trying to create a shapefile with polygons representing areas of different slopes.  I am running into some problems using the following workflow.

I start with a digital elevation model and using Raster>Analysis>Slope I am converting to a slope raster.
I am using "Layer Properties" to assign different colors to different slopes.
So far everything looks right.
Next, I'm using Raster>Conversion>Polygonize.  The result shows many polygons where only one field of color existed in the raster.    I want a single closed polygon for each field of color.

Does anyone know how to achieve this?


texts: How to call for Matrix Parameter / Enum Parameter using Field Calculator in the graphical modeler of QGis?

texts:
I am using QGis Graphical Modeler v.3
It is possible to define, as an Input, a Matrix, defined as follows:
"A table (matrix) parameter for processing algorithms".

I want to use values stored within a Matrix in the Field Calculator, but I do not know how to call for them.

When I open the Field Calculator, I see the reference to other parameters (e.g. Vector Layer, as in the image) but nothing related to the Matrix.

Which is the syntax to call for Matrix values e.g. Matrix (i j)?

Same issue with Enum Parameter (Input).

I searched everywhere but the only related discussion I have found appears to me as related to the GUI: https://github.com/qgis/QGIS/pull/6936


texts:
I previously asked Updating attribute information about this script last week. I have amended it and am confused.

I need to write the ADMINAREA value to the empty LOCALITY fields from TWB_Suburbs to TWB_Property respectively. This script took about 3.5 hours to run and the result did not update the LOCALITY field in the TWB_Property feature class. Can anyone please help?

# Import modules
import arcpy
import os

#Set workspace geodatabase
arcpy.env.workspace = "C:\Users\Harry\Documents\Uni\Semester 1 2016\GIS3407 GIS Programming and Visualisation\Assignment 2\Assignment2.gdb"

#Start an editing session
edit = arcpy.da.Editor("C:\Users\Harry\Documents\Uni\Semester 1 2016\GIS3407 GIS Programming and Visualisation\Assignment 2\Assignment2.gdb")
edit.startEditing(False, False)

# Create a list of suburbs
# Create an empty list
localityList = []

#Create the search cursor
rows = arcpy.SearchCursor("TWB_Suburbs")

#Move to the first row
row = rows.next()

#Append the list with the selected locality name
while row:
    localityList.append(str(row.ADMINAREA))
    #Move to the next row
    row = rows.next()

print(str(len(localityList)) + " <-- # of suburbs")

del row, rows
# Make a layer of all of the properties with no locality value
# Create query
nullQuery = '"LOCALITY" = \' \''

# Create layer
nullLayer = arcpy.MakeFeatureLayer_management("TWB_Property", "TWB_Property_Layer", nullQuery)
field = "LOCALITY"
fc = "TWB_Property_Layer"

#Counter
counter = 0

# Select by location
indexCounter = 0
print localityList
for counter in range(0,len(localityList)):
    localityQuery = '"ADMINAREA" = ' + '\''+ localityList[indexCounter] + '\''
    localityLayer = arcpy.MakeFeatureLayer_management("TWB_Suburbs", "TWB_Sububs_Layer", localityQuery)
    arcpy.SelectLayerByLocation_management("TWB_Property_Layer", "WITHIN", "TWB_Sububs_Layer")

    parcelList = []
    sRows = arcpy.SearchCursor("TWB_Property_Layer")
    sRow = sRows.next()

    while sRow:
            parcelList.append(sRow)
            sRow = sRows.next()
            with arcpy.da.UpdateCursor(fc, field) as cursor:
                for uRow in cursor:
                    if uRow[0] != "str(localityList[indexCounter]":
                        uRow[0] = str(localityList[indexCounter])
                        cursor.updateRow(uRow)
    print ("There are " + str(len(parcelList)) + " empty locality fields in " + localityList[indexCounter])
    arcpy.Delete_management("TWB_Sububs_Layer")
    counter += 1
    indexCounter += 1

del uRow

arcpy.Delete_management("TWB_Property_Layer")
arcpy.Delete_management("TWB_Sububs_Layer")

The output in PyScripter was:

20 <-- # of suburbs

['GLENVALE', 'CRANLEY', 'ROCKVILLE', 'MOUNT KYNOCH', 'HARLAXTON', 'MOUNT LOFTY', 'DRAYTON', 'WILSONTON', 'NORTH TOOWOOMBA', 'NEWTOWN', 'EAST TOOWOOMBA', 'REDWOOD', 'PRINCE HENRY HEIGHTS', 'HARRISTOWN', 'SOUTH TOOWOOMBA', 'CENTENARY HEIGHTS', 'RANGEVILLE', 'KEARNEYS SPRING', 'MIDDLE RIDGE', 'TOOWOOMBA CITY']

There are 1340 empty locality fields in GLENVALE

There are 332 empty locality fields in CRANLEY

There are 1398 empty locality fields in ROCKVILLE

There are 137 empty locality fields in MOUNT KYNOCH

There are 1367 empty locality fields in HARLAXTON

There are 1796 empty locality fields in MOUNT LOFTY

There are 839 empty locality fields in DRAYTON

There are 3608 empty locality fields in WILSONTON

There are 1865 empty locality fields in NORTH TOOWOOMBA

There are 4320 empty locality fields in NEWTOWN

There are 2629 empty locality fields in EAST TOOWOOMBA

There are 109 empty locality fields in REDWOOD

There are 283 empty locality fields in PRINCE HENRY HEIGHTS

There are 3371 empty locality fields in HARRISTOWN

There are 2637 empty locality fields in SOUTH TOOWOOMBA

There are 2448 empty locality fields in CENTENARY HEIGHTS

There are 3421 empty locality fields in RANGEVILLE

There are 2331 empty locality fields in KEARNEYS SPRING

There are 2445 empty locality fields in MIDDLE RIDGE

There are 2379 empty locality fields in TOOWOOMBA CITY


texts:
What are the differences between the various versions (1 through 5) of the Global Land Cover (GLC) 2000 product?

I am going to use northern North America, Europe and Asia (north of 50°N).  Is it better to mosaic the continental products (even though they are different versions) or use the global product (v 1.0).  I want the best data source that I can get for land classification, so if anyone has other recommendations, I would appreciate it.


texts: Using the field calculator to calculate intersections

texts:
I'm comparing several maps of one area (made by different people) to each other, and used intersection to find how much of specific selected features (polygons) overlapped. I soon realised that the intersection did not calculate the overlapping area. I have been told that the easiest way to find the area of the overlap is to do it using the "Field calculator" after the intersection, but I can't figure out how.
If anyone knows an easier way, I'm open to that as well.


texts: How to configure Geonetwork to be accessed from remote machines?

texts:
I just installed Geonetwork and have been playing around with it on our server (localhost:8080), now I want to access this from other machines on the network but when I type in the address I get a page not found error.

I found this link through googling and altered the jetty.xml file (I couldn't find a mapServices.xml file). As a result I can now get to the page but none of the mapping components work (they all still refer to geoserver on localhost). Is there a location where I can make this change once, instead of going through heaps of file and individual changes?

Thanks
Andrew


texts:
I work on a geospatial aggregator. We support certain kinds of temporal data, but we don't have clear names for the different types, and I don't have a clear understanding of the range of geospatial+temporal data. Is there a definitive classification of types? Even a standard term for each one would be helpful.

For instance, these are either supported or planned:

raster + time (eg, satellite imagery)
point in time (eg, earthquakes)
point with time series (eg, river levels at sensor locations)
point moving in time (vehicle locations)

And obviously the above could be applied to polygons and lines, to cover things like change in extent of vegetation, or position of a storm front.

And further, some of the above could be combined, like a point moving in time that also has one or more time series associated with it (eg, a sensor mounted on a vehicle).

So:

Are there other "spatial + temporal" combinations not included in the above?
What are the definitive terms for each combination?


texts: OSM Relation not always in its own Overpass area

texts:
Until now, I assumed that an OSM relation is always in its Overpass area counterpart. For example, the relation 1682152 is returned when searching for administrative boundaries in the area 3601682152, as evident in this Overpass query.
However I have found several relations where this is not the case. For example, the relation 1682200 is not returned when searching in its area 3601682200, as evident in this Overpass query.
Are there specific rules on when an OSM relation is in its own Overpass area (provided the area exists)?


texts: Openlayers 3 - how to keep features from disappear when outside of visible view?

texts:
I have gps coords features on a map and I want to add a line that shows the route of the feature , so I created an icon with long arrow (40px X 600px)
I've noticed that the moment I pan the map and the feature is outside of the visible area , it disappears (although the icon should still be visible)

How can I keep the features from disappear outside the visible map ?


texts: Copying XY to clipboard using Python Add-In (ArcMap)

texts:
I am trying to create a simple tool to copy XY from where i click the mouse button. What i have done already is:

create add-in files using Python Add-In Wizard
create toolbar and tool inside it

in *.py file i got:

import arcpy

import pythonaddins

import win32clipboard as clipboard

class p_tool(object):

def __init__(self):
    self.enabled = True
    self.shape = "NONE"

def onMouseDownMap(self, x, y, button, shift):
    button = 1
    shift  = 2
    clipboard.OpenClipboard()
    clipboard.EmptyClipboard()
    xy = str(x)+' '+str(y)
    clipboard.SetClipboardData(xy, clipboard.CF_TEXT)
    clipboard.closeClipboard()
    message = xy
    pythonaddins.MessageBox(message, "My Coordinates", 0)

About win32clipboard - i tested in in python windon in ArcMap and I am able to import it and openclipboard() but emptyclipboard() does not work and gives me that error:

File "", line 1, in
error: (1418, 'EmptyClipboard', 'Thread does not have a clipboard open.')


texts: qgis composer doesn't lock filter properly

texts:
I'm making a composer with 6 maps based on the same (postgis)layers and extent. The only thing that differs between them is the filter I put on two layers. I change the filter on the layers, then only refresh the map for that particular filter and lock the map. At the end it looks fine in the composer, but after exporting all maps have the same (last applied) filter. This also happens after refreshing the composer.

I'm using a clean install of QGIS 2.18.13 on macOS Sierra 10.12.6


texts: mapserver how to get rid of unknown character from label

texts:
I am developing a gis system using mapserver and open scales library. GIS data is taken from a shape file. When i try to show the name attached to a place its showing unknown characters (not in english) along with its name in english.
here is my .map file settings

 LAYER
   NAME city
   TYPE POINT
   STATUS ON
   LABELITEM 'NAME'
   MAXSCALE 7500000
  DATA "ndia_location.shp"

 CLASS
   EXPRESSION ( '[PLACE]'=='city' )

    LABEL
            COLOR 0 0 0
            TYPE BITMAP
            SIZE SMALL
            POSITION CC
        #MINFEATURESIZE 100
        ENCODING UTF-8
            PARTIALS FALSE
        BUFFER 0
        END
      END

Is there a  way to show only English name for the label.
I have tried with openlayers javascript library also.


texts:
I would like to find local maximum with the moving window 3x3 of a raster image but I don't know how to do this.

Could you please suggest me the way to do this?


texts: Using splitRaster in ArcGIS Pro outputs dark JPEGs, bright PNGs

texts:
I am using the splitRaster function to create 10m x 10m JPEG files from a GeoTIFF, inputting a fishnet grid as a polygon to clip it by. I have figured out the problem of the depth changing between clipping and have set NoData = 254. When I run the tool, however, I am finding that I have very dark .jpg images and very bright .png images as outputs:

^ JPEG example

^ PNG example
I need these images to be in 8-bit JPG or PNG. Ideally, it would look like what is displayed in ArcGIS Pro. I know those images are stretched, but I am unsure how to translate that into these jpgs/pngs.

^ What I would like the images to look like.


texts: ArcGIS Geoprocessing "None Type" Second time

texts:
I am trying to publish an ArcGIS geoprocessing service in 10.2 I have a script setup that will run multiple times through Pycharm, but when I run it in ArcGIS it will run once and return this error:

  Traceback (most recent call last):   File
  "C:\Users\Desktop\git\GeoprocessingService.py", line 477, in
      main()   File "C:\Users\Desktop\git\GeoprocessingService.py", line 372, in main
      cur)   File "C:\Users\Desktop/git\GeoprocessingService.py", line 136, in value_query
      results = cur.fetchall()   File "C:\Python27\ArcGIS10.2\lib\site-packages\psycopg2\extras.py", line
  81, in fetchall
      res = super(DictCursorBase, self).fetchall()   File "C:\Python27\ArcGIS10.2\Lib\decimal.py", line 543, in new
      m = _parser(value.strip()) TypeError: 'NoneType' object is not callable

  Failed to execute (Script).

The script itself is 800 lines, but here is the query I have been running that errors out:

def value_query(id_field, id, value_fields, start_date, end_date, date_field, table, cur):
    query2 = "SELECT "
    for field in value_fields:
        query2 += "NULLIF({0}, '-999')::NUMERIC(6,2),".format(field)
    query2 += "{0} " \
            "FROM {1} " \
            "WHERE {2} = '{3}' AND {0} BETWEEN '{4}' AND '{5}' " \
            "ORDER BY {0} ASC;".format(date_field, table, id_field, id, start_date, end_date)
    #print "This is Query:"
    #print query
    arcpy.AddMessage("{0}".format(query2))
    cur.execute(query2)

    results = cur.fetchall()

    #print "These are value query results:\n\n {0}".format(results)
    return results

Like I said it runs once in ArcGIS then errors out. However, if I close Arc completely and run the exact same script again, it runs once then stops. I feel like there may be some sort of caching or something that Arc does, but am completely stumped.


texts:
I've a shape which was processed with the qgis geoalgorithm tool "eliminate sliver polygons".

My desire is to get a kind of logfile of this certain process. I need to know the whole calculations behind this.

I hope u can imagine what I want..


texts: Can CartoDB provide Web Mapping Service (WMS)?

texts:
Does cartodb implements a wms service?

If not, is there a way to contribute to the project to implement it?


texts: Using IBufferConstruction to create single and multiple buffers (VBA and ArcObjects)

texts:
I need a bit of help with using the iBufferContruction interface, my overall goal is that the user clicks on the map, a graphic is created and then buffered,at multiple distances and then stored in a featureclass

So far, i can create a single buffer (i know if it created because the length that gets returned is correct), however when I store it in my FC the shape is empty.

Here's the code:

Private Sub Buffer_MouseDown(ByVal button As Long, ByVal shift As Long, ByVal x As Long, ByVal y As Long)

'Get the location of the mouse and set it to the point
Set m_pPoint = m_pMxDoc.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(x, y)

'new marker element on the map
Set m_pElement = New MarkerElement

'set the element property
m_pElement.Geometry = m_pPoint

Set m_pGraphics = m_pMap

m_pGraphics.AddElement m_pElement, 0 'Grpahic for the point click.

Set m_pAC = m_pMxDoc.ActiveView

'refresh the map
m_pAC.PartialRefresh esriViewGraphics, Nothing, Nothing
'
Dim pPoly As IPolygon4
Set pPoly = New Polygon

'Dim pPoly As IPolygonElement
'Set pPoly = New PolygonElement

Dim pBC As IBufferConstruction
Set pBC = New BufferConstruction

'Dim pBCProp As IBufferConstructionProperties
'Set pBCProp = pBC
'pBCProp.ExplodeBuffers = True

Set pPoly = pBC.Buffer(m_pElement.Geometry, 200)
Debug.Print pPoly.Length ' RETURNS 1256.63706143592

'pBC.ConstructBuffers
Dim pFeat As IFeature
Dim pfeatureCursor As IFeatureCursor
Dim pfeatureBuffer As IFeatureBuffer

Set pfeatureBuffer = m_pBufferFClass.CreateFeatureBuffer
Set pfeatureCursor = m_pBufferFClass.Insert(True)
Set pFeat = pfeatureBuffer

pfeatureCursor.InsertFeature pfeatureBuffer
pfeatureCursor.Flush
Set pFeat.Shape = pPoly

m_pAC.PartialRefresh esriViewGeography, Nothing, Nothing

Where am I going wrong? Also has anyone got any code examples with creating multiple buffers at various distances?


texts: Leaflet: mutually exclusive overlay layers: "overlayadd" event firing twice

texts:
LeafletJS. I am trying to have two overlay layers be mutually exclusive. Full example here on jsfiddle
I guess it comes down to this bit of code:
        map.on('overlayadd', function(event)
        {
            console.log( "Clicked on " + event.name )
            if ( event.name == "layer1" && map.hasLayer(layer2) ) {
                console.log("remove layer2")
                map.removeLayer(layer2)
            }
            if ( event.name == "layer2" && map.hasLayer(layer1) ) {
                console.log("remove layer1")
                map.removeLayer(layer1)
            }
            ccontrol._update();
        });

It I select just layer2 and then select layer1 it works fine (layer1 becomes selected and shown on map and layer2 toggles to off). The Console shows the sequence of events to be:

Clicked on layer2
Clicked on layer1
remove layer2

But if I select layer1 first, and then try to select layer 2 nothing appears to happen. The console shows the sequence of events to be:

Clicked on layer1 <-- I manually clicked this
Clicked on layer2 <-- I manually clicked this
remove layer1
Clicked on layer1 <-- This happened automatically ???
remove layer2

Can anyone shed some light on this?


texts: Locking/unlocking layout elements using ArcPy with Esri Production Mapping?

texts:
I need to update a text element in a group of MXDs that I manage. That element, however, is locked in the layout using the lock function from the Production Cartography tool.
Using python, is there a way to unlock this text element? In testing, I've found that once it is unlocked I'm able to access the element's content and update it.

Here's the basic code snippet that I'm using:

mxd = arcpy.mapping.MapDocument(md)
for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
    if elm.name == "CompantyTxt":
        elm.text = "New Text to Replace the Old Text"
mxd.save()

I can run this code with the desired text element getting updated if the element is not locked. With it locked, however I get the following runtime error:

  "File "C:\Program Files
  (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects_base.py", line 87, in
  _set     return setattr(self._arc_object, attr_name, cval(val)) RuntimeError: TextElementObject: Error in setting text"


texts:
I'm attempting to debug a 500 server error with my MapServer (running as a FastCGI on IIS 7.5). This would be much less challenging if I could get some indication of where the failure is happening, but MapServer seems disinclined to log anything. I also have no logs in Event Viewer related to the failure so I'm at a loss.

MapServer was configured as documented here. And the map file I'm using is below:

MAP
  NAME           benchmark
  IMAGETYPE      PNG
  EXTENT         -106.39 25.50 -93.31 36.30
  SIZE           400 300
  IMAGECOLOR     255 255 255

  CONFIG "PROJ_LIB" "/MapServer/proj/nad"
  CONFIG "MS_ERRORFILE" "/MapServer/tmp/ms_error.txt"
  DEBUG 3

  WEB
    METADATA
      "ows_enable_request"   "*"
      "wms_title"           "Lonestar Benchmark"
      "wms_onlineresource"  "http://localhost/cgi-bin/mapserv.exe?map=/MapServer/benchmark.map&"
      "wms_srs"             "EPSG:4269 EPSG:4326"
    END
  END

  PROJECTION
     AUTO
  END

  LAYER # States polygon layer begins here
    NAME         ngs_tx
    DATA         "/MapServer/data/tx.shp"
    STATUS       DEFAULT
    TYPE         POLYGON

    CLASS
      NAME       "TX NGS Survey Data"

      STYLE
        COLOR        232 232 232
        OUTLINECOLOR 32 32 32
      END
    END
  END
END

The log directory and file already exist (I first tried without them existing and then created the directory and file after that didn't work). Also, I've tried using "C:" prefixed paths with no difference in behavior.

The CGI is currently set to "connect as" the Administrator, so it can't be a permissions issue. (And yes, this is highly insecure. This is a test environment so doesn't matter.)


texts: Georeferencing metadata in GeoTIFF files

texts:
I am generating maps as TIFFs using the raster package (writeRaster) in R to generate small tiles, then as a final step, combining the tiles via the GDAL script gdal_merge.py.

I thought the resulting files were georeferenced. However, I sent a file to a colleague who is trying to import into ArcMap, and he indicated that the software doesn't recognize the georeferencing.

Here is the result of gdalinfo on one of the problematic files. Can anyone tell me if there's an important metadata component missing, and if so, how best to fix it?

% gdalinfo VERSION18.5_MVN_Vs30_NZGD00_allNZ_AhdiYongWeighted1_noisyT_minDist0.0km_v3_crp1.5.tif
Driver: GTiff/GeoTIFF
Files: VERSION18.5_MVN_Vs30_NZGD00_allNZ_AhdiYongWeighted1_noisyT_minDist0.0km_v3_crp1.5.tif
Size is 11264, 16384
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["GRS 1980(IUGG, 1980)",
        DATUM["unknown",
            SPHEROID["GRS80",6378137,298.257222101],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (1000000.000000000000000,6338400.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 1000000.000, 6338400.000) (166d35'16.17"E, 32d55'40.74"S)
Lower Left  ( 1000000.000, 4700000.000) (165d 1'10.46"E, 47d34'29.83"S)
Upper Right ( 2126400.000, 6338400.000) (178d37'44.88"E, 32d57'57.08"S)
Lower Right ( 2126400.000, 4700000.000) (179d59'25.56"W, 47d38'19.50"S)
Center      ( 1563200.000, 5519200.000) (172d33'57.05"E, 40d28'37.45"S)
Band 1 Block=11264x1 Type=Float64, ColorInterp=Gray
  NoData Value=-1.69999999999999994e+308


texts: Snap! vs. Fulcrum for GIS data collection offline?

texts:
I am evaluating field data collection form-making apps and looking at Snap! or Fulcrum, or something else you know about that I don't!

The key features it needs to have are:

Map-based: I need to load on features to visit and send people out, or alternately let them create a feature as part of filling out a form for it.
Spatially-based records: a user drops a georeferenced point to start a new record and the lat/long/UTM coordinates are part of the attributes.
Required fields: a record should not be allowed to get saved unless required fields have been populated with values.
Pick lists: no-brainer here, both apps do it, but easily managed and updated pick lists are a good thing.
Offline: users need to be offline and still use the map.

...Ideally not a cloud based thing, I am hoping data collectors can export a file (csv, KML, whatnot) which I can open to consolidate multiple data collection dates into a single spreadhseet.  I know both Snap! and Fulcrum depend on cloud syncing, and if I have to go that way so-be-it, but it would be nice to be the owner of my own data.

Any suggestions are welcome, so far my trials on both apps have been so-so.  If any apps out there could take a PDF Form import and turn it into the collection fields/layout that would be even better.


texts: OpenLayers XYZ Unable to Load Local Tiles

texts:
I have some local tiles that were generated from OpenStreetMaps using Maperitive. I've tried going through the docs and tests but the map is not loading. I'm not getting any errors in chrome's js console either.

I think OpenLayers has access to my images because previous errors stated that some tiles could not be found. The solution to that problem was to initialize the OpenLayers.Map.restrictedExtent with bounds. I'm not sure where to go from here.

function initMap() {
    var extent = new OpenLayers.Bounds(-80.604570, 43.228173, -79.481697, 43.607333);
        extent = extent.transform("EPSG:4326", "EPSG:900913");

    var map = new OpenLayers.Map('map',  {restrictedExtent:extent});
    //var options = {'layername':'basic', 'type':'png'};
    var xyz = new OpenLayers.Layer.XYZ('OpenLayersXYZ', 'http://localhost:8081/maps/tiles/${z}/${x}/${y}.png', {sphericalMercator:true, layers:'basic'});
    console.log(xyz);
    map.addLayer(xyz);
    map.zoomToExtent(extent);
}

<!DOCTYPE html>

Maperitive log from tile generation:

set-geo-bounds

Setting the geometry bounds to srid=4326, (-80.604904441331,43.2284165575301,-79.4820319148331,43.6075754238019)

set-print-bounds-view

Set the printing bounds to '-80.604570253079,43.2281730577825,-79.4816977265812,43.6073334440597'

generate-tiles minzoom=0 maxzoom=14


texts: How to repair and clean contours from high data resolution (lidar)?

texts:
I'm unsuccessfully trying to make a bathymetric map with contours extracted from Lidar data at 1 m resolution.
The extracted contours have a lot of artifacts like dangles and islands.

First I used QGIS/GRASS 1.8 with v.clean, running the rmdangle which successfully removed a large part of dangles. But it also created some holes in the main line, and the snapping algorithm don't work well because it's snapping the wrong vertices.

Then, after a few manual operations, the algorithms of Sextante doesn't work anymore and I have an error message : Could not load layer...

So I use ArcGIS 10 for build the topology, but I don't success repairing it. When I try trim or snap, it doesn't work at all.

To conclude I can't remove all the "islands" and other artifacts like intersections.
I think I've tried all I know... And I can't do it manually for the whole study area!

Does anybody have a solution for it?


texts: Clipping raster file using shapefile that has different blocks inside?

texts:
I have a TIFF file as shown below:

I have a shapefile as shown:

How can I separate the TIFF file into different blocks using the above shapefile?


texts: Spatial transformation of latitudes and longitudes that preserve distances between points

texts:
I have a set of lat,lon points.

I would like to apply a spatial transformation (rotate, translate) to place these points somewhere else on the earth surface (for anonymization purposes), but preserve the distance between all pairs of points in the set.

Is there a way to achieve this in Python?

Here is what I have tried, but it gives me errors up to hundreds of meters:

import geopy
import pyproj
from scipy.spatial.transform import Rotation as R

r = R.from_euler('zyx', [10,10,10], degrees=True)

def gps_to_ecef_pyproj(lat, lon, alt):
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    x, y, z = pyproj.transform(lla, ecef, lon, lat, alt, radians=False)
    return x, y, z

def ecef_to_gps(x, y, z):
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    lon, lat, alt  = pyproj.transform(ecef, lla, x, y, z, radians=False)
    return lat, lon, alt

p1 = 51.95103, 12.42487,0
p2 = 51.95103, 13.42487,0
initial_distance = geopy.distance.geodesic(p1,p2)
print('initial distance', initial_distance)

p1_ecef = gps_to_ecef_pyproj(*p1)
p2_ecef = gps_to_ecef_pyproj(*p2)

p1_rot = ecef_to_gps(*r.apply(p1_ecef))
p2_rot = ecef_to_gps(*r.apply(p2_ecef))

final_distance = geopy.distance.geodesic(p1_rot,p2_rot)
print('final distance', final_distance)
print('error', np.abs(initial_distance - final_distance).meters,'m')


texts: Calculating the volume between two interpolation layers in QGIS?

texts:
Given data

From a shapefile with ca. 80 points (displaying borewells) I extracted two interpolations:

interpolation values of the different borewells’s elevations;
interpolation values for the elevation of the bottom of the first
geological layer.

Desired outcome

I need to know the volume of the first geological layer in the area of interest.

I am a beginner and not familiar with GRASS. It would be great if you could explain everything step by step, maybe even providing some screenshots.


texts:
This is my coding I used in Leaflet:

    var pointA = new L.LatLng(28.635308, 77.22496);
    var pointB = new L.LatLng(28.984461, 77.70641);
    var pointList = [pointA, pointB];

    var firstpolyline = new L.Polyline(pointList {
    color: 'red',
    weight: 3,
    opacity: 0.5
    smoothFactor: 1

    });

    map.addLayer(firstpolyline);

The map is not showing any result. I just need straight line from one point to another, just visualize it, not geodetic correct. Points (in Leaflet "circles") are showing perfectly on the map.


texts: NodeJS Geometry string from GPS coordinates before sending to DB?

texts:
I'm building an application that saves geographical data to the postgresql db with postgis extension.

The data I'm saving is in EPSG 4326 WGS84 format, so I'm sending lat/lng coordinates to the postgresql and on insert, it changes the coordinates to a geometry POINT string like so:

GPS (39.10, 9.10)
POINT (1000FF111FFFGGGAAA1) // casual digits and letters

Before sending the data to the db I need to transform the point from EPSG 4326 to 32632 so that I can apply geometric functions to the geometric points. I know how to use ST_Transform, and the translation happens on the DB. What I would like to do is to make the transformation happen on the server so that the point is inserted directly with 32632 SRID.

Is it possible?


texts:
I'm looking for tools and extensions in open source platforms (gvSIG, QGIS, uDig, etc) that allows 3D snapping.

In the last years, most of the projects my company is working on required the editing of 2.5D cartography. 3D snapping is a crucial feature for these projects. And the designers fell more comfortable using CAD tools (like Microstation Geographics) for digitalizing and editing on this projects. ArcGis 9.3 with a 3D snap plugin is used in the second half of the editing part (soon will use ArcGis 10).

Now my company wants to start working with open source platforms (hooray!), and the 3D snapping seems to be one of the main holes in the open source GIS platforms.

Do you know anything could be useful on this subject? Tools, extensions, initiatives, etc.


texts: Creating attribute table for raster data?

texts:
I have raster dataset for rainfall but it does not have an attribute table.

How can I create tabular data using the same raster file and the join back again to it?

Software used is ArcMap 10.1.


texts: ArcGIS - Selecting features depending on 2 attributes

texts:
Ok, I guess this should be not too hard, but I can't figure it out:

I have a layer "A" resulting from a layer "0" which was intersected with itself. So layer "A" contains all intersects.

The relevant columns in the table are:  "Value" and "ShapeArea"  I can see which features(polygons) in "A" are matching with each other as they have the same Shape_Area.
The Intersecting parts have (sometimes) different values in the field "Value" However, what I need are not both intersecting features (polygons)  but only the one with the higher value.

And also I not only have intersections containing 2 polygons but sometimes 3 or more.

As there are 2000 features in the table I can't select them manual...

So in the example in the picture:
Objectid 1 and 2 --> I would only need one of them (which one is not relevant as they have th same value)
Objectid 3 and 4 --> I need Objectid 3 as the value is higher

As a side note I wonder why you can intersect a layer with itself but you can't merge a layer with itself so you could use the combining rules...


texts: Geoprocessing web service not writing out all files for .shp

texts:
I created a Python script, run by a toolbox. I am using Python shapefile to generate a point shapefile. If I run the tool locally (i.e. on the desktop) it works great. Shp, shx and dbf all writen out as expected and the points are valid.

Once I share the geoprocessing script as a service, and then run the tool via an ArcGIS server connection on a remote PC, the only file that writes out is the .shp file. And of course without the .shx and .dbf the .shp file can't be used.

Any ideas on why once the geoprocessing script runs server side, all the files aren't generating for Python library Shapefile?

Any other ideas on how to simply write out a point shapefile using Python?

I know typically a code snippet, is desired, but the code runs locally it's only server side where there is an issue writing the files to a client via, ArcCatalog.


texts: Mapbox data not appearing a lower zoom levels

texts:
I'm styling a Mapbox map, and want to show mountain peak labels. I've created a layer for this, and in the "Select data" tab set the zoom range to 8-20. However on the map (and in the data select preview) the data disappears at zoom levels much higher than 8 - some peaks disappear at zoom 12, some at 10.

Is there some hidden filtering going on and if so how do I control it?


texts: I am looking for Canada Locality/neighborhood data but did not get any useful links. Anybody know about a good Locality/neighborhood data for Canada??


texts: Invalid UUID string in sentinelhub when tryng to access satellite image

texts:
I'm trying to use sentinelhub package on jupyterlab notebook and for some reason I don't manage to be able to acess images.
When I use my credentials in the configuration at the beginning it works and pronts the configuration file location, but whe nI try to acess images as descrribed here I get error message:

DownloadFailedException: Failed to download from: long link of
sentinel-hub    with HTTPError: 404 Client Error: Not Found for url
another long link of sentinelhub
Server response: "Invalid UUID string: my credentials"

This is how I have tried to acess the sat images (copied from the documentation):
!sentinelhub.config --instance_id MySecretCredentials
!sentinelhub.config --show
#that works, print the configurations

from sentinelhub import WebFeatureService, BBox, CRS, DataCollection, SHConfig

INSTANCE_ID = ''  # I have already put instance ID into configuration file at the beginning so left it empty as required
if INSTANCE_ID:
    config = SHConfig()
    config.instance_id = INSTANCE_ID
else:
    config = None

search_bbox = BBox(bbox=[46.16, -16.15, 46.51, -15.58], crs=CRS.WGS84)
search_time_interval = ('2017-12-01T00:00:00', '2017-12-15T23:59:59')

wfs_iterator = WebFeatureService(
    search_bbox,
    search_time_interval,
    data_collection=DataCollection.SENTINEL2_L1C,
    maxcc=1.0,
    config=config
)

for tile_info in wfs_iterator:
    print(tile_info)

what do I do wrong?
My endgoal is to be able to acess sentinel images using sentinelhub on jupyterlab notebook


texts: Spatial Join Attribution for Report Generation

texts:
I generated a grid index for a very large area of land in order to create a mapbook to keep track of some environmental sites for monitoring.  I was able to generate the pages of maps without issue, however, I now wish to tabulate the sites (small polygons) contained within each individual map (grid index cell) in a separate file to append to the map book.  I spatially joined the small polygons to the grid index cells so that I could index the sites by page number.

When I generate a report (extracting only the PageNumber and Site ID in order to create a checklist for each map page), I cannot get the report to group this data into separate pages.  It just generates a running spreadsheet in ascending order.  What I am aiming for is to create a document where every page contains only the sites for each respective PageNumber.
I thought it would be straight forward but I can't figure it out and don't want to export to an excel file to manually separate 75 pages worth of sites.


texts: Align and reshape polygons based on surrounding polylines

texts: I am trying to use a shapefile with polygons and I am trying to align it to a road network, that the polygons align themselves to surrounding streets. Is there a way to automatize alignments between nodes of the polygons and polylines or sometime similar?


texts: Working with spatial data for India when key state truncated due to territory issues?

texts:
I am in India, working with  spatial data. The data is downloaded from gadm.org http://gadm.org/ The state of Jammu and Kashmor (the top most state) is truncated (it reflects ground reality where 1/3 of state is with Pakistan and rest with India).

For legal reasons I need to show the whole state (though obviously I have data for 2/3 provinces, How do I go about it. Almost all global maps show the truncated map which I cannot use.

My tools are R and Python /SAS VA

See below for India:

See below for for GADM image

How can i work (programming) with static map untruncated image and dynamic geospatial (rest of map)?

I have both data sources so its not a data source problem.

How do I merge the two maps (top of one map with bottom of other)?


texts: Suppressing Model-Dialog when running PyQGIS Script

texts:
I wrote a tiny script that has a fixed input and fixed output, i.e. I don't need any user interaction.
However, when running my script the Model-dialog

pops up although the script already ran. Is it possible to suppress the dialog?

I am on QGIS 3.8.3.


texts: Postgis convert long lat from 2 different locations into line and calculate the distance in meters

texts:
The name of my table is odk_duct.duct_monitoring1_core
and the select statement below is what i want to see the distance. The table has 2 longitudes and 2 latitudes.

select
"_URI" as id,
"_SUBMISSION_DATE" as date,
initcap("INFO_TEAM_NAME") as team_name,
initcap("INFO_SUPERVISOR_NAME") as supervisor_name,
"INFO_LOCATION" as location,
"RPT_FTR_POINT1_LAT" as first_point_lat,
"RPT_FTR_POINT1_LNG" as first_point_long,
"RPT_FTR_POINT2_LAT" as second_point_lat,
"RPT_FTR_POINT2_LNG" second_point_long,
"RPT_FTR_WIDTH" as width_cleared,
formid
FROM odk_duct.duct_monitoring1_core

I want to know the distance between point 1 and point 2 in meters using PostGIS query.


texts: ArcMap "Could not open the specified file"

texts:
Has anyone ever encountered this message?  I am getting this when I try to open certain shapefiles.  I am working on a problem that generates shapefiles from spreadsheets using ThinkGeo's .NET tools.  When I try to open ones where I have approximately 142 columns (of mostly string, some float, and a couple datetime fields) I get this message.  I don't get anything else besides this message.

I can open the file in ArcPad 8.0 and QGIS 2.0.1, but it fails every time with ArcMap.  I have run the shapefiles through Andrew Williamson's ShapeChecker, which found no problems.  Going through the Check Geometry tool gives me an ERROR 000229 saying it cannot open the file.  It doesn't have a lock on it, so there is something else happening.

Any tools or suggestions that I should look in to?

Edit:

I am using ArcMap 10.2.2.  I can view the files in ArcCatalog, but when I go to preview them, I get a message saying "Open Failed" as the title and "Error opening feature class".  I tried loading this up in ArcCatalog 9.3.1 and ArcMap 9.3.1.  I am able to see the shapes drawn out, but anytime I try to view the attribute data, it crashes without explanation.

This is point data generated from x/y coordinates that were generated from ArcMap originally.

What I don't fully understand about these errors is how to tell exactly what it is.  I have up to 255 columns of data that I can include.  I am able to include all the data if I break up the columns into smaller chunks.  It seems like for whatever reason my 'no crash' limit is around 142 columns.

As for record count.  ArcMap throws an ERROR 000229 Cannot Open when I try to run it on the shapefile.  ShapeChecker verified that the row count and feature count match though.

I have also checked the field name length as I generate them.  None go over 10 characters.  I haven't found anything definitive yet as to what causes my errors or how to solve them.


texts: I need to do a field mapping and i already have the grid (points). I need to put this points in my phone and use my android device to find this points in the field. I won't have any internet connection.


texts: Split line at regular intervals in shapely

texts:
I have been trying to use Shapely to split a shapefile with multiple features at regular intervals. I have used the split function to do this. I can easily create the points at which the lines should be split, but when I then split nothing happens.

What am I missing?

import fiona
import shapely
from shapely.ops import split
from shapely.geometry import *

schema = {
'geometry': 'MultiLineString',
'properties': {'id': 'int'},
 }

with fiona.open("./test_depth.shp") as inp:
    with fiona.open('result3.shp', 'w', 'ESRI Shapefile', schema) as out:
        for line in inp:
            geom = shape(line['geometry'])
            length = geom.length
            for i, distance in enumerate(range(0, int(length), 10)):
                point = geom.interpolate(distance)
                split_line = shapely.geometry.MultiLineString(split(geom,point))
                print(split_line)
                out.write({'geometry': mapping(split_line),'properties': {'id': 1}})


texts: How can I merge a LOT of polylines automatically?

texts:
I have lots of MultiPolylines, three-dimensional arrays where the bottom level is an array of lat/lngs, the second level is an array of those points (a single polyline), and the top level is the group of lines.  Like this:

[
    [
        [50.12312,0.2541],
        [50.12322,0.2547],
        ...
    ],

        [52.412,1.142141],
        [51.9812,0.21237],
        ...
    ],
    ...
]

The polylines in each group are generally contiguous, and I would like to merge them so that connecting lines within a group become a single line.  How can I do this programmatically without some sort of desktop software?  Is there Python,PHP,JavaScript,etc. code out there for this task or do I have to write it from scratch?


texts: Dynamic Symbolizers w/ CQL not URL encoded

texts:
I am trying user Dynamic Symbolizers in my SLD:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd">
<NamedLayer>
    <Name>cite:test</Name>
    <UserStyle>
      <Title>Test</Title>
        <FeatureTypeStyle>
          <Rule>
            <PointSymbolizer>
              <Graphic>
                <ExternalGraphic>
                  <OnlineResource xlink:type="simple" xlink:href="http://localhost/icons/${some_type}"/>
                    <Format>image/png</Format>
                </ExternalGraphic>
              <Size>35</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

How if the value of 'some_type' has characters in it that need to be URL encoded the URL after parsing the CQL is invalid.

An easy example is if the value has a space in it i.e. 'Some Value'.  Which produces the URL http://localhost/icons/Some Value, and it should probably produce: http://localhost/icons/Some%20Value.

Am I using CQL correctly in this instance.  I looked through the provided functions here:

https://github.com/geotools/geotools/tree/master/modules/library/main/src/main/java/org/geotools/filter/function

and did not find one specific to URL encoding.  That being said its not a huge deal to write a filter function to do this.  However if you can use CQL in the href of an OnlineResource the core CQL parsing logic should probably handle this.

Is there anything else I might be missing?  Crossing my fingers ;)


texts: How can I hide only the left panel in Google My Maps?


texts: Exporting an image to Firebase Storage/Google Drive from Firebase Cloud Function

texts:
I could export an image to Google Drive and Cloud Storage from Earth Engine code editor. But exporting image to Cloud Storage from server side doesn't work when included in a Firebase Cloud Function. (earthengine node client library is installed)
Note: Server Side Authentication with a Service Account is successfully done.Followed all the steps mentioned in NPM installation of GEE


texts: Using errorsalm() from spdep with multilevel data?

texts:
I have a dataset that has multiple observations per county, and I am interested in estimating the effect of a county-level variable on individual-level outcomes. I would like to account for adjacency between counties in this model. I am using errorsarlm() from package spdep. It seems errorsarlm() will only accept data frames that have the same number of rows as the length of the listw object.  Since the listw object is the length of the number of unique counties, and the data frame has multiple observations per county, I receive the following error:

Spatial autoregressive error model
Error in errorsarlm(ss_form, data = temp, listw = w, method = "eigen",  :
  Input data and neighbourhood list have different dimensions

Can someone please help me figure out a way to fit this model?

I have provided example code below.

#Get spatial polygons data frame from rgeos
spdf <- getData('GADM',country='USA',level=2)

#Create unique county IDs that will match SSA ID's
spdf$State_County <- paste(toupper(spdf$NAME_1),toupper(spdf$NAME_2),sep='')

#Create neighborhood object
nb <- poly2nb(spdf, row.names = "State_County")

#Create listw object
w <- nb2listw(nb, style = "B")

#Estimate model with multiple observations per county (this is just a mockup of the model I am fitting, it will not actually run)
model1 <- errorsarlm(injury_outcome ~ county_ambulance_use + age , data=trauma, listw=w, method="eigen", quiet=FALSE)


texts: Importing multiple layer files like NAS/OSM in a standalone PyQgis script

texts:
I am trying to convert a dlm50 NAS file of Rheinlandpfalz (Germany) to shapefile. With ogr2ogr this repeadetly converts Polygons to Lines and whole layers are missing. If I import it into QGIS and use the validate algorithm it works fine though (ogr convert format still fails).

Importing them in QGIS gives a pop-up where I can chose all layers I want to import. I don't find a way to do this in a python script though. Yes I can open them with QgsVectorLayer and ogr driver but I receive one object, which can't be used properly.
I had also problems with this with OSM-files and spatialite files which have failed to load with the spatialite driver but succeded with ogr import (may have been to big?).

I know I should try to convert the dlm50 file to PostGIS but it will be the same driver for reading and I am generally interested if I miss something and the import of such files is possible somehow in a standalone python script (not in QGIS Python console). It should also be possible to work with such files without the need to preliminary convert them to other formats?!

So how can I import layers of multiple layer files like NAS/OSM/Spatialite in a standalone python script?


texts: Plotting Points in R with Google Map context. Error in download.file(urlStr, destfile, mode = "wb", quiet = TRUE)

texts:
I am trying to plot point data with Google Map context but i keep getting an Error. I tried the same using OSM context and it worked just fine.

R Code below

library(RgoogleMaps)
Lat <- as.vector(quakes$lat)
Long <- as.vector(quakes$long)
MyMap <- MapBackground(lat=Lat, lon=Long, zoom = 10)
PlotOnStaticMap(MyMap,Lat,Long,cex=tmp+0.3,pch=1,
                col= '#FB6A4A80')

Error Report

"https://maps.googleapis.com/maps/api/staticmap?center=-24.655,176.9&zoom=4&size=640x640&maptype=mobile&format=png32&sensor=true&key="
Error in download.file(urlStr, destfile, mode = "wb", quiet = TRUE) :
  cannot open URL 'https://maps.googleapis.com/maps/api/staticmap?center=-24.655,176.9&zoom=4&size=640x640&maptype=mobile&format=png32&sensor=true&key='
In addition: Warning message:
In download.file(urlStr, destfile, mode = "wb", quiet = TRUE) :
  cannot open URL 'https://maps.googleapis.com/maps/api/staticmap?center=-24.655,176.9&zoom=4&size=640x640&maptype=mobile&format=png32&sensor=true&key=': HTTP status was '403 Forbidden'


texts:
I am using "ST_Intersection" between two 3D lines which I believe do not intersect, but the result gives me their intersection.

SELECT ST_ASText(ST_Intersection(
        ST_GeomFromText('LINESTRING( 1 1 1, 5 5 1)'),
        ST_GeomFromText('LINESTRING( 4 1 10, 1 5 10)'
      )));

results in:

  POINT Z (2.71428571428571 2.71428571428571 5.5)


texts: Geoprocessing tools fail though CRS same, prj-file ok, edit mode off, topology valid, etc

texts:
For my layers none of the QGIS geoproccesing tools (difference, intersection, etc. all give empty attribute tables) work, though I've checked for all (?) error sources (as in the title). Then I tried to do a manual edit of my polygon-geometry (adding a ring) and this is what happend (the data is in my Dropbox):

I'm using QGIS 2.10.0 on Win7


texts: Python/PostGIS - Type "Geometry" Does Not Exist

texts:
I am a total novice with both PostGIS and Python - I am just trying to get this game running. I have followed the instructions here and here and am now at a dead end.
I have pgAdmin 4 downloaded and running with PostGRESQL10 and the PostGIS extensions. However, when I try to set up the database with Python in Command Prompt (file location>python tools\set_up_db.py), I receive the following message:

Connecting to database
        ->host='localhost' port='5432' dbname='sessions' user='postgres' password='xxxxx'
CREATE TABLE IF NOT EXISTS dggrid (id BIGSERIAL PRIMARY KEY, gid bigint UNIQUE, geo geometry, center geometry, population int, employment int);
Traceback (most recent call last):
  File "tools\set_up_db.py", line 22, in <module>
    cursor.execute(query)
psycopg2.ProgrammingError: type "geometry" does not exist
LINE 1: ...(id BIGSERIAL PRIMARY KEY, gid bigint UNIQUE, geo geometry, ...
                                                             ^

Looking around, it seems like this may be an issue with my Schema. My default Schema is "public," and I've tried to change the database's search path to "public, postgis," which hasn't helped. No idea if I'm on the right track here at all. As far as I can tell, I've done everything else correctly in the instructions...


texts: Calculating Moran's I on large point dataframe in R

texts: I have a dataset of about 77,000 points of which I want to test autocorrelation. I have tried to compute the test statistic using the testSpatialAutocorrelation tool in the DHARMa package, which uses "a simulation-based approach to create readily interpretable scaled (quantile) residuals for fitted (generalized) linear mixed models." However when I try to run the test, I run into a memory error. What is a way to easily calculate Moran's I for a large dataset?


texts: How to select multiple segments along one roadway?

texts:
I have a roadway network shapefile which only contains the center line of the roadway. I'm trying to create directional roadways from the center line. So I use "Copy Parallel" to create roadways at both sides of the center line and then add attribute to specify the direction. This is not a good method since the length of the lines do not change, so the ends of each segment can't meet (See 1st pic). Is there anyway to do it like tracing?

Another problem is that after adding the parallel lines, I need to specify the direction of the parallel lines (whether it is northbound or southbound roadways). Since one roadway is consist of multiple segments, it takes a lot of time to manually select roadway of one direction. I can't use "Select By Attribute" since it is difficult to differentiate the parallel lines after copying and pasting.So I'm wondering if there is a way to select the whole road (multiple segments) by picking the start and end point of the roadway (So all segments directly connect to each other will be selected)?

Any suggestions would be very much appreciated.


texts:
I'm trying to get to a piece of code that lets me print PDF maps with the PYQGIS API using a workspace and a template.

This code does the job, but throws an error, stopping me from lining up consecutive maps to be produced.

sys.path.append('/Applications/QGIS.app/Contents/Resources/python')
import os
from qgis.core import (
    QgsProject, QgsComposition, QgsApplication, QgsProviderRegistry)
from qgis.gui import QgsMapCanvas, QgsLayerTreeMapCanvasBridge
from PyQt4.QtCore import QFileInfo
from PyQt4.QtXml import QDomDocument
from qgis.core import *

def make_pdf():
    project_path = 'project.qgs'
    template_path = 'template.qpt'
    gui_flag = True
    QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/Resources/python", True)
    app = QgsApplication(sys.argv, gui_flag)
    app.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)
    app.initQgis()
    canvas = QgsMapCanvas()
    QgsProject.instance().read(QFileInfo(project_path))
    bridge = QgsLayerTreeMapCanvasBridge(
    QgsProject.instance().layerTreeRoot(), canvas)
    bridge.setCanvasLayers()
    template_file = file(template_path)
    template_content = template_file.read()
    template_file.close()
    document = QDomDocument()
    document.setContent(template_content)
    composition = QgsComposition(canvas.mapSettings())
    composition.loadFromTemplate(document, {})
    composition.refreshItems()
    composition.exportAsPDF('output.pdf')
    QgsProject.instance().clear()

The output is:

QSqlDatabase: QSQLITE driver not loaded
QSqlDatabase: available drivers:
QSqlQuery::prepare: database not open
ERROR: Opening of authentication db FAILED
WARNING: Auth db query exec() FAILED
Loading a file that was saved with an older version of qgis (saved in  2.12.1-Lyon, loaded in 2.16.3). Problems may occur.
QObject::startTimer: QTimer can only be used with threads started with    QThread

Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

That SIGSERV error at end seems pretty nasty - but I'm no where near understanding what's going on when I Google that.


texts: SAGA 2.2.5 Problem Object Based Image Segmentation

texts: Every time I use the Object Based Image Segmentation with Landsat 8 data I get a segmentated image which looks like a spider's web. The image contains no 0 values. Anyone has a idea how to fix this?


texts: Monitoring layer loading progress in OpenLayers

texts:
I'm doing an upgrade from OpenLayers 2 to OpenLayers 3. In OL2, I was monitoring the loading progress of layers with their loadstart and loadend events. With grid layers I was even able to show a progress bar using tileloaded event.

In OL3, I can't find such events, neither in documentation, nor in the code.

Do you know about some way to at least catch the moment when loading or rendering of a layer is finished?


texts: STContains, STIntersects and STWithin return wrong result for geography loaded from WKT


texts: ZoomTo not working in ArcGIS Runtime SDK for WPF?

texts:
I have a graphics layer. When the user selects an address from our software I want the map to ZoomTo that address. This works, but the user has to have an interaction for the map first. I thought this might be because the Map Extent was being reset, but when stepping though the code, it seems to Initialize first. Below is the code I am using. I separated the zooming code into another method just to see if it made a difference, but it did not. I also tried taking the extent off of "MyMap"

public void ShowMap(string account_num, string street, string city, string state, string zip)
        {
            Locator FindlocatorTask = new Locator("http://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer");
            FindlocatorTask.AddressToLocationsCompleted += ShowMapLocatorTask_AccountToLocatonCompleted;
            FindlocatorTask.Failed += ShowMapQueryTask_Failed;
            AddressToLocationsParameters addressParams = new AddressToLocationsParameters();
            Dictionary<string, string> address = addressParams.Address;
            string formatAddress = " ";
            formatAddress = street + " " + city + ", " + state + " " + zip;

            OStreet = street;

            if (street != null)
            {
                RelatedFeatureGrid.Visibility = Visibility.Visible;
                OptionsPortal.Visibility = Visibility.Collapsed;
            }

            accountAddress = account_num; //used to find the associated meter
            RelatedAddressTextBox.Text = formatAddress;

            if (!string.IsNullOrEmpty(formatAddress))
                address.Add("Address", formatAddress);

            FindlocatorTask.AddressToLocationsAsync(addressParams);

        }

        private void ShowMapLocatorTask_AccountToLocatonCompleted(object sender, ESRI.ArcGIS.Client.Tasks.AddressToLocationsEventArgs args)
        {

            CandidateListBox.Items.Clear();
            _candidateGraphicsLayer.ClearGraphics();

            try
            {
                AddressCandidate candidate = args.Results[0];
                Graphic graphic = new Graphic();
                graphic.Symbol = LayoutRoot.Resources["DefaultMarkerSymbol"] as ESRI.ArcGIS.Client.Symbols.Symbol;
                graphic.Geometry = _mercator.FromGeographic(new MapPoint(candidate.Location.X, candidate.Location.Y));
                graphic.Attributes.Add("Address", candidate.Address);
                string latlon = String.Format("{0}, {1}", candidate.Location.X, candidate.Location.Y);
                graphic.Attributes.Add("LatLon", latlon);
                _candidateGraphicsLayer.Graphics.Add(graphic);

                ZoomToLocation(graphic);
            }
            catch
            {
                MessageBox.Show("Address Cannot be Found", "Error", MessageBoxButton.OK);
            }
        }

        public void ZoomToLocation(Graphic graphic)
        {

            MapPoint candidatePoint = graphic.Geometry as MapPoint;
            double displaySize = MyMap.MinimumResolution * 30;

            ESRI.ArcGIS.Client.Geometry.Envelope displayExtent = new ESRI.ArcGIS.Client.Geometry.Envelope(
                    candidatePoint.X - (displaySize / 2),
                    candidatePoint.Y - (displaySize / 2),
                    candidatePoint.X + (displaySize / 2),
                    candidatePoint.Y + (displaySize / 2));

            MyMap.ZoomTo(displayExtent);
        }

No errors are being thrown.

I actually think I know what it is but I have no idea how I am going to handle it. It seems that the Initial Extent of MyMap is being set AFTER the Extent of the ZoomTo. The ZoomTo basically beats the initialize event, even though it is fired first.

I have tried wrapping the code in an "if (MyMap.IsInitialized)", but it still is a no go. HOWEVER, if I put a messagebox in the ZoomTo method that forces the user to click ok, it works every time. I really don't want to have to do that.

Thread.Sleep also does not work because it stops the Map from loading completely


texts: Road names of intersections in a city to CSV - Overpass API

texts:
I would like to get the road names of intersections in a particular city in Overpass API.
I already have some code to get the nodes of the intersections, but I could not solve to get the road names which meets. I would like the output to be CSV file.
http://overpass-turbo.eu/s/10oo code here or here:
//[out:csv ("id","name1","name2","traffic light")][timeout:250];
{{geocodeArea:Budapest}}->.searchArea;
way[highway~"^(motorway|trunk|primary|secondary|tertiary|(motorway|trunk|primary|secondary)_link)$"](area.searchArea)->.major;
way[highway~"^(unclassified|residential|living_street|service)$"](area.searchArea)->.minor;
node(w.major)(w.minor);
out;


texts: Geoserver : create a rasterized polygon / dynamic fill for polygons with features from database

texts:
I am using geoserver 2.7-SNAPSHOT, everything is working fine but I need to implement a new function or tranformation process for personnal needs. I have a map with differents base layers (provinces,neighboorhood, district...etc) and an other layer which appears upon the base layers. This layer use a dataset which is a score, corresponding to the influence of a store on a certain area. The score is from 0 to 1.07. I succeed to create a choropleth map

The problem is that I got one color per district (each polygon is a district) so I tried to create an heatmap but then the division isn't visible, I have a heatmap for the whole country.

What I'd like to do is having either an heatmap or at least a dynamicfill with gradient depending of the value of the feature inside each district, not a static color. So there will be different gradient to each district

To do so I tried to use a new function, using the source code for heatmap, colorramp, dynamicfillcolor that I found on geotools repositorybut I failed.

How can I implement a new function to my geoserver ?

I succeed to add a new process (and it was hard since the documentation related to geoserver is poor), but it doesn't work for a function.


texts: Reorganizing OpenStreetMap into PostGIS database entity?

texts:
How do I reorganize the OSM given entity to more flexible and schematic purpose?

The OSM gave 4 entities the planet_osm_line, planet_osm_point, planet_osm_polygon and planet_osm_roads. Is there a good practice to reorganize them  like into city_building, city_road and etc. I'll be reorganizing a city's spatial objects.


texts: ArcGIS Javascript Map loads but doesn't display

texts:
I just started creating a very simple mobile GIS app using ArcGIS Javascript API, and I have experiencing a very peculiar issue.
When the app loads, the basemap is white (not visible) until you use the zoom control.  After that, the app works except for the BasemapToggle control.  I was using another similar post as a template to add a custom basemap to the BasemapToggle control seen here:
How to Use Custom Basemaps with BasemapToggle widget?
There are no errors reporting, so I am not sure what is causing both the fault basemap load and the BasemapToggle issue.
Has anyone else experienced this and have information on resolving the issue?
Any help appreciated!


texts: Passing multiple Raster layers from C# to ModelBuilder?

texts:
I have created a model in ModelBuilder and I need to pass a list of rasters from my code (C#) but until now, ModelBuilder gets only the first record in the array.

I have been passing it in different ways like:

IVariantArray parameters = new VarArrayClass();
parameters.Add(lstRoots);//lstRoots is a list of strings with the root of the raster.

Other way:

foreach (string r in lstRoots) //lstRoots is a list of strings with the root of the raster.
{
   parameters.Add(r);
}

And another way:

IVariantArray parameters = new VarArrayClass();
IVariantArray paramlist = new VarArrayClass();

foreach (string r in lstRoots)
{
   paramList.Add(r);

}
parameters.Add(paramList);
result = (IGeoProcessorResult)GP.Execute("ModelCellStatistics", parameters, null);

In the model I have created a variable that is the parameter with multiple values:


texts: QGIS 3.10 problems with using the "Reshape features" tool

texts:
Theoretically, I know how to use the Reshape tool but in practice, I don't know what I am doing wrong.

I am double-clicking on the node, upon it's become active (green color), then creating some new nodes, until the finish. Afterward, I am doing the right-click and everything is gone, or the shape doesn't change the borders as I want when I do the loop selection.

Could someone explain to me how to use this tool properly?

Below the situation:

My proposal new border

And everything is gone after the right-click, the same as after escape button.

In turn, when I do something like the loop,

My shape is shrinking, whilst I want to have it extended (towards up according to the image).

UPDATE:

I can use the Vertex tool feature, but it's still far from my goal, as I have to modify all separate nodes one by one.


texts: ArcSDE 10.0 giomgr Error During Post Install

texts:
I just installed ArcSDE 10.0 and ran the post install once to test the install and it ran great.  Today I tried to run another post install to create another database and I get the error:

"Parameter "TEMP" not found in giomgr.defs file!"

This occurs in the repository setup.  Has anyone had this same error or something similar? A point in the right direction would be great!


texts: How to Get Client Credentials When User is Already Logged Into Server Directory (ArcGIS JavaScript 3.25)?

texts:
I'm working on an ESRI JavaScript App.  My app needs to support OAuth2 authentication.

The code below shows how the client logs in and how I get the client's credentials.  What I need is the token variable from the client credentials.

The code below works when the user isn't logged in.  The credentials array comes back as empty when the client is logged into the ArcGIS Server endpoint.  What's really strange is that the code works when the client is logged into our Portal.  It just fails when the client has logged into our ArcGIS Server Directory.  The reason I need the token is that I'm sending GET requests as queries to certain feature layers.  The get requests have to have a token themselves since the services are secured.
How can I get the client credentials when the client is already logged into ArcGIS Server?

    //create the OAuth2 handle for logins
    var info = new OAuthInfo({
        appId: "appID",
        portalUrl: "portalURL",
        popup: false
    });
    esriId.registerOAuthInfos([info]);

    //get token for signed in user
    esriId2=esriId.toJson();
    token=esriId2.credentials[0].token;


texts: Is srsName of "EPSG:4326" spec compliant in GML 3.1 / 3.2?

texts:
I sometimes see srsName="EPSG:4326" in GML documents (e.g. https://postgis.net/docs/ST_GeomFromGML.html, although my particular source is not public). The intent (except perhaps for axis order that isn't important here) is clear, and probably interoperable across various systems, but I wonder if that is actually specification compliant.

I've looked at the OGC Naming Authority and OGC URN policy, and they use the urn:ogc:def:crs:EPSG:\d{4,5} form (i.e. it would be urn:ogc:def:crs:EPSG:4326).

The GML Point Profile (OGC 05-029r4) uses srsName="urn:ogc:def:crs:EPSG:6.6:4326" in examples, and requires that to be "well known". I don't see anything in there that actually excludes any particular srsName though.

Is the srsName="EPSG:abcd" form (in particular, srsName="EPSG:4326") specification compliant for a GML 3.1 or 3.2 element (in particular, for Point Profile)?


texts: How to evaluate the statistic significance in the Neighbour Nearest Test?

texts:
I'm reviewing a scientific article to be publish.
I tried to evidence that pathogens collected from different geographic areas clustered. I used Nearest Neighbour test to prove it. However, I've noticed that no p-value is assigned in my output. So, how could I get the p-value from my analysis? It will be for a publication. Otherwise, how could I explain the significance of the result?

This is the output obtained from my analysis using QGIS3:

Observed mean distance: 495.1212145321877

Expected mean distance: 654.1975182291224

Nearest neighbour index: 0.7568375005036007

Number of points: 28

Z-Score: -2.46153773279898

Could anyone help me?
Thank you in advance


texts: Update image geometry from another geometry (GEE)

texts:
I have a set of images I have processed. Unfortunately, the geometry of the image processed do not always fit the image data.

Problem: Let's assume I have an imageCollection. Then, when I would like to filter by bounds according the green area (using filterBounds). The 'blue' image will be selected whereas I don't want (no data in the green area).

Illustration:

blue is the image data
dark the image geometry
green is the region of interest

Question: Knowing the geometry where there is data (I know it from image before pre-processing), how can I refit the image geometry ?

My attempts:

to update the property system:footprint with:

image = image.set('system:footprint', mygeometry)

But, then, when I call image.geometry(), the geometry hasn't been updated.
clip function: the documentation clearly explains this is applying a mask, do not update the geometry.

Reproducibility: Link code


texts:
I have two coordinates A = (60◦ N, 30◦ E) and B = (60◦ N, 150◦ E). I want to calculate the shortest possible route between these two points and also I want to know how far north does the route goes.

I have no idea how to calculated the how far north does the route goes so I am asking help.


texts: Too many features in QGIS DXF import: 1 Feature per Column

texts:
The company I work for do a considerable amount of field surveys, and export the data directly from the GPS (Leica GS08) to DXF files.
Whenever I try to open a DXF file (with lines and points), the points in QGIS are displayed in a strange fashion.

I opened the same files in ArcGIS, and the result was distinct:

As you can see, the total amount of features is quite different. In ArcGIS, the system loads 48 features (the same number of points present on the GPS survey), while QGIS loads 288 individual features.

In fact, as it is visible on the Attribute Table, QGIS loads 6 features per GPS point: one feature per "Coordinate", "Height", "Attribute", "Point_Code", "Attribute" and "0" (which are records of the "Layer" column), while ArcGIS loads all these attributes, in one feature.

I have imported the file by dragging the DXF into the QGIS system. I also tried "Another DXF Import/DXF2Shape" plugin, but the result was not the desired one. Although the image presented only one point, there were not a database associated, and the label of the text was not practical.

My issue is that I need QGIS to open these dxf files on the same way ArcGIS does - where each feature may have distinct columns on the attribute table, but not a feature per attribute. As one can imagine, if there is a big concentration of data, is almost impossible to work with the file.

Please find attached a DXF example of the situation in hand. DXF File Example


texts: Python: Trying to figure out Geopy in python for reverse look ups

texts:
I am trying to test the following code on my system but I cannot get the correct output.

Is this code correct for using geopy with GoogleV3?
If so, how can I return country and continent name?

import subprocess
from geopy.geocoders import GoogleV3

def set_proxy(proxy_addr):
    cmd = "set HTTPS_PROXY=" + proxy_addr
    call = subprocess.call(cmd, shell=True)
    cmd = "set HTTP_PROXY=" + proxy_addr
    call = subprocess.call(cmd, shell=True)

set_proxy("my.proxy")

api_key = "my key"
geolocator = GoogleV3(api_key)
location = geolocator.reverse("52.509669, 13.376294", timeout = 10)
#would like to return just country and continent name


texts: How to add temporary PostGIS table to canvas?

texts:
I'm working with PostGIS database. I would like to create temporary table with geometry and add it as a layer to QGIS. There is no problem with regular tables, but I can't figure out how to add temporary tables. I checked source code of DB Manager plugin and found this key line: uri.setDataSource("", u"(%s\n)" % sql, geomCol, "", uniqueCol). With this in mind I created my code:

uri = QgsDataSourceURI()
uri.setConnection(server, port, db_name, login, password)
sql = "CREATE TEMPORARY TABLE temp_1 AS SELECT * FROM some_spatial_table"
uri.setDataSource('', u'(%s\n)' % sql, 'geom', '', 'gid') # 'geom' and 'gid' are present in the table
vlayer = QgsVectorLayer(uri.uri(), 'layer_name', "postgres")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

The SQL-query is Ok - I'm able to execute it in DB Manager and add resulting table to map. But my code doesn't produce a layer in QGIS and there is no exception raised. What is wrong?

UPD:
There is an output in PostGIS tab in QGIS logs:

Unable to execute the query.
            The error message from the database was:
            ERROR:  syntax error at or near "CREATE"
            LINE 2:                         CREATE or replace TEMP VIEW sp_temp_...
                                            ^
            .
            SQL: SELECT * FROM (
                                    CREATE or replace TEMP VIEW sp_temp_2015_10_25_02_29 AS
                                    WITH posts AS (
                                    SELECT post_name, date, pollutant, measurement, type_measure
                                    FROM air_observations
                                    WHERE pollutant = '301' AND date = '2014-05-05'
                                    )
                                    SELECT air_posts.gid, air_posts.geom, air_posts.name,
                                    posts.date, posts.pollutant, posts.measurement, posts.type_measure
                                    FROM air_posts
                                    INNER JOIN posts
                                    ON air_posts.name = posts.post_name;
            ) AS "subQuery_0" LIMIT 1

Starting from "CREATE" and ending with ";" is an actual query that is passed as sql.


texts:
I am trying to iterate through some feature classes, named "Hos_AllPatients2014", and select the certain data and name it "Hos_IP", then export it to the corresponding folder "Hos". I have already made folders for each "Hos" that contain another layer I had to make for each. However, there are some "AllPatients2014" feature classes for which I do not have a folder. I thought I deleted them all out, but found out halfway through my model running that it had an error and stopped. I would have liked to have had the model contain a portion to check if "Hos_IP" already exists, and if so move on to the next iteration. Instead, I had to create a new folder with the remaining "AllPatients2014" and started the model again.

I would like to find out how to get the model to make this check for a feature class and skip to the next iteration if it exists. I would also be interested to figure out if the folder for "Hos" exists, and if not, create a folder for it and go ahead and export the data like the rest.

I don't know how helpful a picture of my model will be, but it is included. Also open to hearing how I can refine the whole thing.


texts: Do I need PostGIS to get Points in range?

texts:
I'm going to build a table of places where each contains geographical coordinates (latitude and longitude).

In my basic request to the database, I want to get all records in range X km to the location with specific lat and lng.

For example, get the restaurants in the range 50 km from point: { lat: 12.345, lng: 67.891 }.

Is it right idea to use PostGIS for this kind of query? Or can I use plain SELECT query without PostGIS features?


texts: Getting area (in square meter or km) of polygon in MySQL 5.6


texts: I installed a fresh version of the QGIS LTR 3.10. I am working on creating a postgres server from a raspberry pi4 and connecting the database to my qgis. Through a tutorial I am walking through, they suggest creating a Basic Authentication configuration (image 1). I type in a username and password, click the convert to configuration button and I receive a box asking for the current master password (image 2). I've tried researching how to find this and have only come across opening the qgis-auth.db file. I've opened it with a text editor, but its mainly a bunch of symbols that don't make sense to me. So I am wondering how/where do I set up this master password?


texts: Geometric Network - Listening to network feature connection changes in the IObjectClassEvents::OnChange event

texts:
I need to be able to tell when a particular type of feature, that participates in a geometric network, is having its geometry automatically modified due to the geometric network.

What i mean is that I'm currently listening to the IObjectClassEvents::OnChange and use it to work out if the feature being changed is of the type I am interested in and if the change is to its geometry. I use IFeatureChanges and IRowChanges. This works fine. However I also want to check if the change is caused by a change to a connected node in the geometric network OR if the change is due to the user modifying the feature directly.

Anyone have any idea how to do this ?

I would be looking for a INetworkFeatureChanges or an INetwokClassEvents interface, but they don't exist..

Thanks,
Ian


texts: lidR:::catalog_laxindex() - "largest interval gap increased to..." logged

texts:
Often when running lidR:::catalog_laxindex() on a catalog of tiles I'm seeing the following log messages:

largest interval gap increased to 1554
largest interval gap increased to 3802
largest interval gap increased to 15080
largest interval gap increased to 3921

while the function is running. The function runs successfully to completion, so should I worry about these messages or ignore them?


texts: ArcMap not enforcing one to one relationship classes

texts:
I'm designing a data survey application.

I created a fire hydrant feature class named devBorne_geom. This feature class contains only the geometry and cannot contain survey fields as I want them to be independent from the main feature class

I added a devBorne_attrib table containing all the attributes i need for my app. It does contain a BorneID field to be linked against the devBorne_geom.GlobalID field

I defined a relationship class between these two tables. Its defined as composite, one-to-one, with Forward notifications.

My problem is that, in arcMap(happens in arcGIS Pro too), when I'm editing devBorne_geom features, I'm able to add multiple devBorne_attrib records for a given devBorne_geom feature, despite having a one-to-one relationship.

Is there something I'm doing wrong?

EDIT

Following J.R's suggestion, in ArcCatalog, I customized the relation like this, enforcing a 0 to 1 relation between the origin and the destination

This modification didn't prevent from adding multiple records, but ArcMap's Validate Feature tool did flagged features having more than one devBorne_attrib. I now have to figure out how to do it in ArcGIS Pro


texts: What projections should I use to make my own Globe?

texts:
Searching for an answer to this question, I found  instructions posted by Gulf of Maine Research Institute showing how to create a globe.

Using manual methods ...

What approach would I take to create a globe using GIS?

What projection should I use for each individual gore?

If I wanted fewer seams near the poles, is there some other projection I could use?

Could I do a series of projections to create the gores based on a soccer ball and stitch them together?

How would I determine the point of tangency for each pentagon and hexagon, along with their vertices in latitude/longitude?

Would some other non-soccer-ball isohedron be more suitable?


texts:
I am trying to map some planetary data on my webMap server using OpenLayers. I defined a proj4 projection as follow:

Proj4js.defs["MYDATUM"] = "+proj=longlat +a=1740000 +b=1740000 +units=m +no_defs";

In the html file I also call the proj4js.js, and then I added the following lines:

 map = new OpenLayers.Map('map', {
       controls: [
               new OpenLayers.Control.Navigation(),
               new OpenLayers.Control.PanZoomBar(),
               new OpenLayers.Control.ScaleLine({bottomOutUnits: ''}),
               new OpenLayers.Control.MousePosition(),
               new OpenLayers.Control.Graticule()

               ],
       Projection:  new OpenLayers.Projection("MYDATUM")
  });

My issue is that the scaleline does not show the correct distance and seems not to take into account my own projection. Actually I tried a couple of example and it seems that the line:

   projection:  new OpenLayers.Projection("MYDATUM")

does not have any effect on the ScaleLine size. How can I display correctly a scale bar using my own projection/datum ? Again, I'm not working with terrestrial data, so I do not want wgs84 or any of it's friends.