texts: ArcMap Crashes with ArcGIS has encountered serious error when Accessing Classified Symbology

texts:
Whenever I select classified symbology for a raster layer ArcMap will crash/close and give the message/ prompt to send error report: "ArcGIS has encountered a serious error.."

The problem is not associated with the data itself as it happens for any raster.

I have uninstalled ArcGIS 10.3 and reinstalled ArcGIS 10.4 and reset my normal file and the error still happens. There doesn't seem to be any issues with the install and the issue is new. I had been using ArcGIS 10.3 for a while with no issues.

The only change on my computer I can think of is I updated my AMD Radeon Display driver recently. However I have tried disabling it and running ArcMap and the same error occurs.


texts: Making polygons slighty bigger QGIS 3.14.16-Pi

texts:
I am not much of a programmer, but I have been using QGIS for about a year now to do fibre planning. I use the Atlas print to PDF feature to print the suburb plans with a block overlay layer to print individual blocks of a suburb. Around all the blocks there lines representing fibre lines, and on these lines there are points representing access points. In this week I started using  rule based symbology and labelling where I would use the following code in each layer to turn off all layers except the layer that is active on the atlas:
within($geometry, @atlas_geometry )=1
This works great if the polygon covers everything in the block in question, including the lines and points around the block.
My blocks where generated by a Model in QGIS that was written by someone else, and only covers the erven within the block, and not the surrounding lines and points, which I need to print as well. So this is what I am getting in my print layout at the moment. Only the addresses and the block's label:

This is what I need:

All Addresses are showing, including the erven lines, the lines around the block and the point around the block with their labels.
The problem is that all 128 blocks I need to export to PDF as individual pages, are all too small. I need to slightly enlarge each polygon to include the points and lines around it.
I have tried affine transform, buffer vectors and transforms vector layer without success. I read that the grass v.transform might be able to perform this function, but it gives an error when I run it:


texts: Merging raster data with different resolutions

texts:
I need to merge different raster datasets which have one layer and contain discrete values to a layer stack.
The resolutions of the data can have any value so not only multiples of the resolutions from the other input rasters thus there are problems with alignment. Also the corner coordinates are not necessarily identical.

Has anybody got a creative idea of what is the best way to do the merging without losing information as a layer stack can only have layers with the same resolution? I have used gdal_merge.py and of course there is quite some loss of information through the nearest neighbour resampling (though I must use nearest neighbour as I want to keep discrete values).

One thing I thought of was to polygonize each raster and then union them and rasterize the resuling polygon with a appropriate resolution. That would keep most information but here I am afraid that it will not work properly with large datasets. Perhaps somebody has experience with this?


texts: NLCD WMS Server Returns Color Index Instead of Classification Code

texts:
When I use the WMS service to download a region of the NLCD 2016 dataset the service returns a geotif with color index value and not classification code.  The image successfully loads in ArcMap or QGIS and is rendered using the appropriate style.  However, when I identify individual cells to check their values I get color index and not classification code.
Furthermore,  I can use the online viewer and download a similar image using the draw tool.  When I do this the downloaded image contains the correct style and classification codes instead of color index. This is desired output I am looking for from WMS service.
Should my query be adjusted to return the correct value (classification) instead of color code?
Ex. query -
https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2016_Land_Cover_L48/ows?version=1.3.0&service=WMS&layers=NLCD_2016_Land_Cover_L48&styles&crs=CRS:84&format=image/geotiff8&request=GetMap&width=226&height=100&BBOX=-101.462828,35.636649,-101.456214,35.638981&


texts: How to select and add a featureclass through GXDialogue

texts:
I am new to arcobject
i have learned how to add a feature class by specifying its name in code through gxdialogue.

i want to change it so i could directly select the the featureclass instead mentioning its name in coding (i.e. want to remove hardcoding).

can you please suggest what changes i need to do in code.

below is the code:

  IGxDialog gxd = new GxDialogClass();
        gxd.AllowMultiSelect = false;
        gxd.ButtonCaption = "Add featureclass";
        gxd.Title = "Add FeatureClasses";
        gxd.RememberLocation = true;
        IGxObjectFilter gxObjFilter = new GxFilterPersonalGeodatabases();

        gxd.ObjectFilter = gxObjFilter;
        IEnumGxObject gxEnumObj;
        gxd.DoModalOpen(ArcMap.Application.hWnd, out gxEnumObj);
        IGxObject gxObj = gxEnumObj.Next();
        //getting the address of PGDB
        string WSAddress = gxObj.FullName;

        IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;

        IWorkspaceFactory wsf = new AccessWorkspaceFactoryClass();

        IWorkspace ws = wsf.OpenFromFile(WSAddress, arcMap.Application.hWnd);

        IFeatureWorkspace fws = ws as IFeatureWorkspace;
        IFeatureClass fc = fws.OpenFeatureClass("FetureClassName");

        IFeatureLayer fl = new FeatureLayerClass();
        fl.Name = "Feature";
        fl.FeatureClass = fc;

        mxdoc.AddLayer(fl);
        mxdoc.ActiveView.Refresh();
        mxdoc.UpdateContents();


texts: Adding Route Shapefile (PolylineM) to QGIS?

texts:
I'm new to QGIS (and to this amazing forum as well) and I already have an issue :

I received a Route Layer in SHP format (PolylineM). I've easily loaded this file into QGIS but it doesn't seem to recognize the "M" aspect of it.

I couldn't find anything recent about it but I've read that in the past, QGIS couldn't deal with PolylineM. Is it still the case?

Problem is I have nothing to re-calibrate that Route if that's the only option.


texts: How do I load mapinfo mid/mif into geopandas?

texts:
I suspect this is a trivial question to answer but how do I load mapinfo mid/mif files into geopandas?

While it's trivial to load and manipulate .shp files, I can't work out what code to use:

Here's a few examples of the data I'm trying to play with: http://www.stoke.gov.uk/ccm/content/business/general/open-geospatial-consortium-data-catalogue.en

cylpaths = geopandas.GeoDataFrame.from_file(path)

I tried to test the geometry, which shows:

cylpaths.crs

 Out[15]:
{u'ellps': u'airy',
 u'k': 0.9996012717,
u'lat_0': 49,
u'lon_0': -2,
u'no_defs': True,
u'proj': u'tmerc',
u'towgs84': u'375,-111,431,-0,-0,-0,0',
u'units': u'm',
u'x_0': 400000,
u'y_0': -100000}

But when I try to preview the data, I get an empty dataframe:

cylpaths

Out[16]:
Int64Index([], dtype='int64')   Empty GeoDataFrame

I'm a bit lost in working out what I need to do next.


texts: Checking for the existence of a workspace (arcpy), issues with if else statement

texts:
I am trying to describe, list, and test for existence of featureclasses in a geodatabase using the following code.

import arcpy

#Set geoprocessing environments
arcpy.env.workspace = "C:""

#Check if the environment workspace exists
if arcpy.Exists(arcpy.env.workspace):
    print "Workspace exists" + " " + "..." + " " + "getting coordinate systems" + " " + "..."
#If workspace does not exist, statement will print "does not exist"
else:
    print "Workspace does not exist"

fcList = arcpy.ListFeatureClasses("*", "POLYGON")

for fc in fcList:
        desc = arcpy.Describe(fc)
        print fc + "" + ":" + "" + str(desc.SpatialReference.name)

While the current code does work when the workspace does exist, when I try to give it a non-existing workspace I get an error message:

The issue appears to be around the line code that allows me to print the list of feature classes in the geodatabase along with there spatial reference information
for fc in fcList:
TypeError: "NoneType" object is not iterable

I simply want it to print statement following else: to print "the workspace does not exist", but I do not know how to actually accomplish that.


texts:
When reading data from a VRT, gdal appears to interpolate values from all available rasters. Is there a way to configure the VRT in such a way that it selects the first valid raster in the dataset, or is there a different raster format that would allow this?
I'm looking to migrate from ossim to gdalwarp, and ossim uses an ordered list of datasets when looking up elevation data. gdal's RPC_DEM doesn't appear to support multiple dataset sources like ossim, so I've created VRTs for all our elevation datasets and one VRT combining all of them. Gdal, however, seems to interpolate between all available rasters, including low quality data that should only be used as a last resort, introducing unnecessary error:
$ for f in $(ls *.vrt); do echo $f; gdallocationinfo -wgs84 -valonly $f -117.5 33.75; done
1.vrt
458
2.vrt
456
3.vrt
458
4.vrt
460
5.vrt

6.vrt
621
merged.vrt
605

I've looked into the documentation and experimented with SimpleSource but didn't get anywhere. I think I can emulate it by masking each VRT in the list with the inverse of the previous VRTs, but that seems like it would be awfully complex, inefficient, and it may not work correctly near mask boundaries.


texts: ArcGIS symbology mis-categorizing data points

texts:
Wondering if anyone has run into this issue with Arc before. I've set up a simple symbology categorized into groups by parcel size.

Size Categories:

This appears to be calculating correctly for the most part, however a handful of data points are being styled as "Other" when they should be placed in a specified category. For example this data point is being styled as if it were out of range when parcel size attribute is within 0-10 acres.

I'm mostly confused as to why this is an occasional thing, if none of my points were rendering correctly I would assume that there was a problem with my data or with the symbology classes I set, however like I mentioned earlier this is only happening for some tens of points out of thousands.


texts: Reformatting API Call results in an Atlas. Nominatim and OpenWeatherMap

texts:
I am creating information sheets for geotagged photos using the atlas features in QGIS 3.14. Some of the information is created using a graphical model but I want to add some more using API calls.
I have managed to use HTML frames to present the results of calls to Nominatim and OpenWeatherMap using expressions like the below
'https://nominatim.openstreetmap.org/reverse?format=xml&lat='||  "latitude"  || '&lon=' || "longitude"  || '&zoom=18&addressdetails=1'

This works well but the results contain more information than I need so I would like to be able to re-format them.
<reversegeocode timestamp="Wed, 07 Apr 21 09:25:16 +0000" attribution="Data © OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright" querystring="format=xml&lat=88.88&lon=88.88&zoom=18&addressdetails=2">
<result place_id="101970880" osm_type="way" osm_id="46739254" ref="Gypsy Lane" lat="51.6895477895927" lon="-0.43032246110111033" boundingbox="51.6858356,51.69023,-0.4312647,-0.4303209" place_rank="26" address_rank="26">Gypsy Lane, Hunton Bridge, Abbots Langley, Leavesden, Three Rivers, Hertfordshire, East of England, England, WD4 8PS, United Kingdom</result>
<addressparts>
<road>Gypsy Lane</road>
<hamlet>Hunton Bridge</hamlet>
<suburb>Abbots Langley</suburb>
<village>Leavesden</village>
<city>Three Rivers</city>
<county>Hertfordshire</county>
<state_district>East of England</state_district>
<state>England</state>
<postcode>WD4 8PS</postcode>
<country>United Kingdom</country>
<country_code>gb</country_code>
</addressparts>
</reversegeocode>"

The only strategy I can think of is to try and download an xml for each photo in the graphical model then read the contents into a string and re-format it but that seems clumsy and I probably couldn't get it to work.


texts: Can QGIS read/load a zipped shapefile or Geopackage from a URL/URI?

texts: Does anyone know if it is possible to load a vector layer in QGIS (not stored in a RDBMS) from a URL/URI?


texts: Postgis calculating distances between polypoints on big dataset

texts:
I have a table that looks like this:

CREATE TABLE public.gps_data
(
  gpstimestamp timestamp with time zone,
  geom geometry(Geometry,27700),
  elevation double precision,
  accuracy double precision,
  bearing double precision,
  speedms double precision,
  subjectid character(100),
  id integer,
  curveid integer,
  distance double precision,
  usable integer
)

This is filled with about 6 million poly points.
My aim is to calculate distances between those points, whereby the distance stored per row will be the distance of the current point to the point that comes before it if ordered by gpstimestamp or id.

Right now I'm doing that with the following query from php:

SELECT ST_Distance(a.geom, b.geom) FROM testdata a, testdata b WHERE     a.id=" . $i . " AND b.id=" . ($i - 1)

With $i running through the number of rows I have.
Using this method I have to do 1 query per row which of course takes forever (on my laptop i get about 2000 queries per minute).


texts: Is it possible to update attribute from spatial join without creating new shapefile?

texts:
I have the following shapefiles:

county.shp: polygons referring to country land subdivisions
lake.shp: polygons referring to lakes, ponds, etc.

Layers attributes are as follow:

county.shp

id INTEGER
name TEXT
code INTEGER

lake.shp

id INTEGER
code INTEGER

The codeattribute is not filled for lake.shplayer. I would like to update this field by intersecting both layers and bring back codevalues from county.shp to lake.shp.

I know I can do this with the Join Attributes by Location tool, but this produces a new shapefile with columns put together. Is it possible to do a simply update lake.shp layer based on a spatial intersection?


texts: How to read csv and compare with a shp to create a point layer with these attributes?

texts:
I need to read a csv and compare with a polygon shapefile. If some cell combine with a feature attribute of the shp, This feature is added to a point layer, with a correspondent row with data of the csv and the centroid of the polygon shp. In the code below I try to do that but didn't work. How can I do this?

from qgis.PyQt.QtGui import *
from qgis.PyQt.QtCore import *
import processing

source_layer = iface.addVectorLayer("/path/polygon.shp", "Setor", "ogr")
crs = source_layer.crs().toWkt()

# create a point layer
centroid_layer = QgsVectorLayer('Point?field=cd_geocodi:string&crs=' + crs, 'MyPoint' , 'memory')

# read a .csv
uri = "path/teste.csv?delimiter=%s" % (",")
lyr = QgsVectorLayer(uri, 'DataBase','delimitedtext')
lyr.fields()[0].name()
QgsMapLayerRegistry.instance().addMapLayer(lyr)

centroid_layer.startEditing()

for feat in lyr.getFeatures():
    for source_feature in source_layer.getFeatures():
        geometry = source_feature.geometry()
        centroid = geometry.centroid()
# select the attribute
        name = source_feature.attribute("CD_GEOCODI")
        if feat[0] == name:

# Create the new feature with the fields of the memory layer
# And set geometry and attribute before adding it to the target layer
            print feat[0], "match with", name
            centroid_feature = QgsFeature(centroid_layer.fields())
            centroid_feature.setGeometry(centroid)
            centroid_feature['cd_geocodi'] = name
            centroid_layer.addFeature(centroid_feature)
        else:
            print feat[0]

        centroid_layer.commitChanges()

# Add the layer to the registry
QgsMapLayerRegistry.instance().addMapLayer(centroid_layer)


texts: Unable to load and display data from geoJSON source in OpenLayers

texts:
It seems like OpenLayers library changes from day to day and loses compatibility with the previous minor version each time... But, actually, I have another question to ask. The problem I face, is that I'm unable to load data from url. In previous versions I did it like so:

var vectorSource = new ol.source.GeoJSON({
  projection: 'EPSG:3857',
  url: '../assets/data/countries.geojson'
});
var vectorLayer = new ol.layer.vector({
  source: vectorSource
 });

Now, there are no longer any ol.source.GeoJSON class and no ol.layer.vector class (with ol.layer.Vector instead). So, how to implement it? The second problem I face, is that I'm unable to create and display vector layer from local geojson object. I know how they do it now and I did it in the same manner:

var geojsonObject = {
"type": "FeatureCollection",
"features":[{
    "type": "Feature",
    "geometry": {
        "type": "Point",
        "coordinates": [82.0, 0.5]
    },
    "properties": {"prop0": "value0"}
},{
    "type": "Feature",
    "geometry":{
        "type": "LineString",
        "coordinates": [
            [102.0, 0.0],
            [103.0, 1.0],
            [104.0, 0.0],
            [105.0, 1.0]
        ]
    },
    "properties": {
        "prop0": "value0",
        "prop1": 0.0
    }
},{
    "type": "Feature",
    "geometry": {
        "type": "Polygon",
        "coordinates": [
          [ [100.0, 0.0],
            [101.0, 0.0],
            [101.0, 1.0],
            [100.0, 1.0],
            [100.0, 0.0] ]
        ]
    },
    "properties": {
       "prop0": "value0",
       "prop1": {"this": "that"}
    }
}]
};

var vectorSource = new ol.source.Vector({
  features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
});

var vectorLayer = new ol.layer.Vector({
  source: vectorSource
});

map.addLayer(vectorLayer);

But when I run this code, I get only some default blue circle in the very centre of the map (no errors and no warnings in the console, but still some wrong behaviour).

I wonder what is the right way to do it?


texts: In what data format should I export a (lat,long,alt) or (x,y,z) grid to be compatible with the largest number of GIS software?

texts:
I am new to GIS and developing an application that generates a 3D grid with (latitude, longitude, elevation) or (x, y, z) metric triplets, and I would like to add an exporting feature that ensures compatibility with the most GIS software.

In case of several possibilities, I'd love to understand the pros and cons considering the current industry standards and also the "efficiency" of the format, as I'd like to reach the largest number of users.

UPDATE: I did a bit of homework but I am still a bit lost with the variety of formats. What I have come accross which seem interesting: GRD (Surfer), KML (Google), DXF (AutoCad), WRL (seems cool but not used anymore), SHP (ArcView - closed format?)

Anyone could help me choose?


texts: ArcGIS JS API: Basic map not loading inside dojo BorderContainer

texts:
I am facing this issue while loading the simple basic arcgis map inside dojo BorderContainer. Following is the code that I am trying,

this.mapContainer = new BorderContainer({
    region: 'center',
    design: 'sidebar',
    style: 'height:calc(100%-150px);border:0px solid #b7b7b7;width:100%;overflow: hidden;'
});
this.mapContainer.startup();
this.mapContainer.layout();
this.mapContainer.placeAt(this);

dcMap = new Map(this.mapContainer.domNode,{
    basemap: "topo",
    center: [-122.45, 37.75], // longitude, latitude
    zoom: 13
});

But, when i load the app i get following outcome

So not only some tiles are missing from the map, the zoom control as well not working properly.

At first i tought the issue was with dom i.e. the container might not be ready when the map tries itself to load in it. But it seems like thats not the case, cause if I put the timeout and load map after some delay I get the same results.


texts:
I'm using now for a couple of years Augmented Reality to visualize 3D-Models, but I never thought to much about the correct projection, when I create my 3D-Models. But for the current project I'm working on I noticed some gaps between 3D-Models which should actually be together.

Currently I'm using a self programmed QGis-Plugin to convert 2D-LineStrings into 3D-Pipes. The first step is to convert the 2D-LineStrings into a metrical coordinate system, in my case a UTM coordinate system.
After that I do my calculation for the vertices, faces and normals to get a 3D-Pipe. Then I save every single 3D-Pipe to a .obj/.mtl file.

As an AR-SDK I'm using Wikitude, but I noticed the same behaviour as described above with Metaio (which is not available any more).

I take my 3D-Model and place them with a WGS84 coordinate into the AR-Scene.

The anchor point is normaly the first point from the 2D-LineString.

Roughly it's working but not exactly. Probably I get better results if I'm using the pseudo mercator coordinate system, but I'm not sure.

I'm also not 100% sure what kind of procetion is used in an AR-Scene. I guess that it is a metrical system (OpenGL) which is maybe transformed to a global system like WGS84, which would explain why I could place 3D-Models with a WGS84 coordinate.

To recap everything:

In which metrical coordinate system should I transform my 2D-LineStrings to do my calculations to get the most accurate results in an AR-Scene.
What kind of projection is used in an AR-Scene.


texts: Difference between mapped unit and target unit

texts:
I am looking to use the SoilGrids REST API to fetch soil properties.
The coordinate inputs for lon-lat are (-72,-9). Some of the soil properties show two units of measurement - mapped units and target units. Can anyone please help me in understanding why I am seeing two different units for one property? Which one is correct and what is the rationale behind showing two units?
Sharing a snippet from the API response below.
{
  "name": "ocd",
  "unit_measure": {
    "d_factor": 1000,
    "mapped_units": "g/dm³",
    "target_units": "kg/dm³",
    "uncertainty_unit": ""
  },


texts:
I add a shapefile of the district in Paris on another shapefile. I would like to match them. See images below:

The black border is the original, the blue one is what I want to get.

I tried with the editing polygon function, espacially reshape polygons, offset x,yand then the delete parts function. The last one deleted what I wanted :

But now I need to close them. I could maybe draw them but I would like to rely on the shapefile below for more accuracy.

So :

How to make new border based on the shapefile underneath?
Maybe there is another way?
How could I draw them if there is no other way ?

Using @snaileater's answer worked as I wanted , I had  to reshape the original shapefile (first image) with the vertex tool  and to choose in the snapping options current layer.

It was a bit laborious because it was manual, so I was asking myself if there would be a way to do it with a mask or a buffer ?


texts: Difference between EPSG:4326 and EPSG:32637?

texts:
I'm exporting data from a program into a .csv for RStudio, which will then go into ArcMap.
I haven't done this before and I've been given two projection options which are:

WGS 84 [EPSG:4326]
WGS 84 / UTM zone 37N [EPSG:32637]

What is the difference between the two of them and how do I know which one I need to use?


texts:
I tried to connect to MSSQL from QGIS and I couldn't find any solution to my problem.

I install MSSQL on my Mac with Docker and try to connect QGIS with localhost(host=localhost) and left the DSN empty, also inserted USER&PASS.

I tried to insert the IP instead, even mentioning the port (1433), nothing helped.

Did anyone have an idea how to solve this issue?


texts: Duplicating feature layer with ArcObjects?

texts:
After opening a MxD document and getting the main map then a feature layer, I want to duplicate the layer with ArcObjects in order to inherit its properties (symbology, layer def and labelling). Then, I want to edit the name of the copy and its layer defs query.

I managed to get both IMap and source feature layer.

IMap map = mapDocument.getMap(0);
IFeatureLayer sourceFeatureLayer = map.getLayer(0);
IFeatureLayer destinationFeatureLayer = ...?
destinationFeatureLayer.setName("layer_2");
map.addLayer(destinationFeatureLayer);

How do I duplicate the source layer?


texts: Examples and materials about creating real time web gis application?

texts:
Nowadays all I see when I explore more about web gis applications are some examples of interactive maps which are actually not so 'dynamic' as they are called. It's usually just simple interactive map with few base layers and overlays with markers and popups on click. If they are little better then they offer you some kind of filtering and querying the data.

It's good way to start learning how to develop web gis applications but what about real time web gis applications?

It looks like very interesting field and more and more possible with new frameworks and technologies.
E.g. let's say that you want to have map which shows you movement of dolphins which are tagged with satellite chip and that chip sends location data to web server (e.g. each 20sec) and then it's display on map in real-time? That's the first thing that come to my mind but there are hundreds of possible examples.
I'm really struggling to find good examples of such projects.


texts: bar charts on MapServer are not plotted depending on the zoom level

texts:
I have a problem with bar charts on mapserver. At first glance the charts are perfectly plotted and all is okay. But depending on the zoom level a few charts are not plotted which are plotted on other zoom levels.

The problem appeared on MapServer Version 5.6.5. Afterwards I tried Version 6.4.1 and the problem remains (but other charts are vanishing). Also the geometry type of the source (polygon or point) makes no difference.

Is there something wrong with the Mapfile or should I not use Shapefiles as source? Here is the Layer-Definitition:

LAYER
   NAME "Gebiete_p"
   GROUP "Gewerbegebiete_p"
   TYPE POLYGON
   STATUS ON
   DATA GewGeb/GewGeb_Gewerbe_Privat.shp
    CLASS
        NAME "Gem"

        STYLE
           COLOR 120 120 120
           OUTLINECOLOR 20 20 20
           ANTIALIAS FALSE
           OUTLINEWIDTH 1
           #WIDTH 1
        END # STYLE
    END # Class

     PROJECTION
            "init=epsg:25832"
     END

         METADATA
                "wms_title" "Gewerbe_p"
                "wms_srs"   "EPSG:25832 EPSG:32632 EPSG:31466 EPSG:4326 EPSG:31468 EPSG:31467 EPSG:4314 EPSG:4258 EPSG:900913"
                "wms_server_version" "1.1.1"
                "wms_format" "image/png"
         END

    END #LAYER

    LAYER
       NAME diagr3_p
       Group Gewerbegebiete_p
       TYPE CHART
       #DATA GewGeb/GewGeb_Gewerbe_Privat.shp
       DATA GewGeb/GewGeb_Gewerbe_Privat_Zentroid.shp
       STATUS ON
       #LABELCACHE ON
       PROCESSING "CHART_TYPE=bar"
       PROCESSING "CHART_SIZE=35 45"
       #PROCESSING "CHART_BAR_MAXVAL=70"
       #MAXSCALE 14999
       #MAXSCALE 44999
       #MINSCALE 35000
       #SYMBOLSCALEDENOM 1000000

       CLASS
         NAME "1 Mbit"
         STYLE
           SIZE [p_diag1]
           COLOR [p_farbe1]
           OUTLINECOLOR 0 0 0
         END
       END
       CLASS
         NAME "2 Mbit"
         STYLE
           SIZE [p_diag2]
           COLOR [p_farbe2]
           OUTLINECOLOR 0 0 0
         END
       END
       CLASS
         NAME "6 Mbit"
         STYLE
           SIZE [p_diag6]
           COLOR [p_farbe6]
           OUTLINECOLOR 0 0 0
         END
       END
       CLASS
         NAME "16 Mbit"
         STYLE
           SIZE [p_diag16]
           COLOR [p_farbe16]
           OUTLINECOLOR 0 0 0
         END
       END
       CLASS
         NAME "30 Mbit"
         STYLE
           SIZE [p_diag30]
           COLOR [p_farbe30]
           OUTLINECOLOR 0 0 0
         END
       END
       CLASS
         NAME "50 Mbit"
         STYLE
           SIZE [p_diag50]
           COLOR 0 0 255 #[p_farbe50]
           OUTLINECOLOR 0 0 0
         END
       END

     PROJECTION
            "init=epsg:25832"
     END

         METADATA
                "wms_title" "diagr3_p"
                "wms_srs"   "EPSG:25832 EPSG:32632 EPSG:31466 EPSG:4326 EPSG:31468 EPSG:31467 EPSG:4314 EPSG:4258 EPSG:900913"
                "wms_server_version" "1.1.1"
                "wms_format" "image/png"
         END
    END

On the first zoom level the large grey area have no charts ...

on the second it have ...

on the third the chart in the middle is suddenly vanished


texts:
I have a map image that I have reprojected from OSGB36 to googleMaps projection (EPSG:900913). When I run gdal2tiles on the reprojected map, the results are not alligned with the google maps images (its not just out by a little bit but by 10's of miles).

If I run it through MapTiler (from maptiler.org) it is aligned fine. MapTiler asks for an SRS which I dont specify when using gdal2tiles. So I'm thinking I've missied out a stage?

Also when I run gdalinfo on the warped file, the coords of the corners seem wrong. E.g. for the bottom-left corner the latitude is 51.143906, but it should be 50.955887 for google maps.

I'm usign gdal16 installed by OSGeo4W on Windows7 64-bit.

Warped file here..
https://docs.google.com/leaf?id=0B93lksnTC7_cNzIzZWQ0MDYtMTY5Yy00NjY0LThiMjktZjc2YzA3MmU5NjI4

gdalinfo on warped file...

C:\OSGeo4W\apps\gdal-16\bin>gdalinfo e:buckland_WARPED.tiff
Driver: GTiff/GeoTIFF
Files: e:buckland_WARPED.tiff
Size is 611, 419
Coordinate System is:
PROJCS["Google Maps Global Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235630016,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-474051.292871501410000,6617009.987366309400000)
Pixel Size = (8.392224865731517,-8.392224865731517)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -474051.293, 6617009.987) (  4d15'30.51"W, 51d 9'49.57"N)
Lower Left  ( -474051.293, 6613493.645) (  4d15'30.51"W, 51d 8'38.06"N)
Upper Right ( -468923.643, 6617009.987) (  4d12'44.69"W, 51d 9'49.57"N)
Lower Right ( -468923.643, 6613493.645) (  4d12'44.69"W, 51d 8'38.06"N)
Center      ( -471487.468, 6615251.816) (  4d14'7.60"W, 51d 9'13.82"N)
Band 1 Block=611x3 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=611x3 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=611x3 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=611x3 Type=Byte, ColorInterp=Alpha

The sequence of commands I run are...

Georeference the OSGB PNG images...

gdal_translate -a_srs EPSG:27700
-gcp 0 400 241538.5 119812
-gcp 600 0 244713.5 121928.66666667
-gcp 0 0 241538.5 121928.66666667
-gcp 600 400 244713.5 11981
bucklands1.png bucklands1.tiff

Warp into googleMaps projection (EPSG:900913) ...

gdalwarp -s_srs EPSG:27700 -t_srs EPSG:900913 -of GTiff
bucklands1.tiff bucklands1_warped.tiff

Tile (make into googleMap 256x256 tiles) ...

gdal2tiles bucklands1_warped.tiff gmTiles

It works if I by-pass the gdalwarp stage and make google tiles directly from the OSGB36 images. But the image quality is poor as the interpolation is not so good.


texts: How to catch PostgreSQL/PostGIS errors in GDAL/OGR Python bindings?

texts:
Is it possible to catch PostgreSQL output error messages in OGR after a query has failed?

For example

conn.ExecuteSQL('SELECT * FROM non_existing_table;')

In Postgres this returns an error message:

"Error: relation 'non_existing_table' does not exist."

In ogr we get the standard:

"Error: current transaction is aborted, commands ignored until end of transaction block"

It would be a bonus if we could catch/parse the Postgres error message. I suspect this is a bit of a long-shot given what I've read about SWIG/CPL, but worth asking to see if anyone's got ideas. Ideally this would be in Python, but I'll take C/C++ comments as well!

Thanks.


texts: Using spectral library provided by USGS Spectroscopy Lab?

texts:
In the USGS Spectroscopy laboratory, the researchers there have worked very hard to measure reflectance of a hundreds of minerals (and plants).

Look at this site

I downloaded the splib04a - FTP Compressed Binary but I don't know how to access its entries in MATLAB.

How can one read and access every spectrum in this spectral library?


texts:
I want to display a map from MapServer in Openlayers 2.13. I've prepared the Mapfile:

MAP
    NAME "WMS Server"
    IMAGETYPE PNG
    EXTENT 106.637 -6.36199 106.78 -6.231
    SIZE 600 300
    IMAGECOLOR 255 255 255
    SHAPEPATH "../data"
    FONTSET "../fonts/fonts.txt"
    UNITS dd
    WEB
        TEMPLATE 'example_template.html'
        IMAGEPATH "C:\ms4w/tmp/ms_tmp/"
        IMAGEURL "/ms_tmp/"
        MINSCALEDENOM   1000
        MAXSCALEDENOM   250000
        METADATA
            "wms_title" "WMS Server Kota Tangerang Selatan"
            "wms_onlineresource" "http://localhost:8000/cgi-bin/mapserv.exe?mode=map&map=c:/ms4w/apps/coba1/test3/wms1.map&"
            "wms_enable_request" "*"
            "wms_srs" "EPSG:4326"
        END
    END
    PROJECTION
        "init=epsg:4326"
    END
    LAYER
        NAME "jalan1"
        METADATA
            "wms_title" "Jalan1"
        END
        # DATA jalan1
        TYPE    LINE
        STATUS  ON
        CONNECTIONTYPE POSTGIS
        CONNECTION "host=localhost dbname=db_mapserver user=admin_map password=admin port=5433"
        DATA "geom from jalan1"
        MINSCALE 1000
        MAXSCALE 100000
        PROJECTION
            "init=epsg:4326"
        END
        CLASS
            NAME 'jalan1'
            STYLE
            COLOR 255 0 0
        END
    END
END

I tried to access the WMS from QGIS it showed the layer, but when I added it to the map it shows nothing, like this picture:

My assumption is there must be something wrong with the Coordinate System, though when I tried to add the data directly from Postgis it displayed correctly, please help me.

regards


texts: How to create a raster in GDAL with Int64 datatype: Rasterizing a vector layer with very large ID as attribute

texts:
I use the following python code to rasterize a vector layer using gdal.RasterizeLayer

output = 'F:\\TEMP\\rasterized.tif'
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Int32)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))

gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=id"])

The id used as the ATTRIBUTE is larger than Int32.

I retrive the layer from a database with:

conn = ogr.Open('DB Connection parameters')
sql = "SELECT Statement"
layer = conn.ExecuteSQL(sql)

Notes:

The GDAL reference does not mention an Int64 datatype, which is consistent with the GeoTIFF format description.
However, the format description also mentions, that BigTIFF is supported when running libtiff >= 4 (which is my case). Reading the libtiff website I could not find a specific reference if Int64 is supported or not.
I installed gdal 2.3.2 via miniconda

Questions: How to create a raster in GDAL that supports the Int64 datatype? Does not have to be GeoTIFF.


texts: Getting cell value along Flow Direction using ArcPy?

texts:
I have:

a flow direction raster where I took off (assigned NoData value) the cells that overlayed the river network (from flow accumulation).
a perfectly overlapping raster where I asigned a coefficient value to the cells (except the ones overlapping the river network)

I need to calculate an index for every cell of my raster that is calculated as the mean between the coefficient value of the cell and those of all the cells between the first and the river network, along the flow direction.

To explain better:
Let's begin with a cell, who has its coefficient value (c1) and has a 8 in the flow direction raster. I want to pick the coefficient value and "move" to the south-west cell (that is the next one along the flow direction). Then I want to pick its coefficient value (c2), see in the FlowDirection raster which is the next cell and move there, to pick its coefficient value (c3)....and so on till I find a cell with NoValue (river network).

Then I need to calculate the mean beetween these values and assign the result, in a new raster, to the beginning cell

(c1+c2+....cn)/n


texts: What are the differences between a WFS and a simple API?

texts:
So I have worked with OSM in the past and have been really happy with the API of Overpass turbo. Now I want to integrate some other data and that one is provided via WFS.

Why is some geodata offered in WFS and other via API?

Is there any difference for me as a Python developer between the two services?


texts: polytiles.py not showing large areas of water with osm carto style

texts:
I've used polytiles.py with the old XML OSM style to generate tiles from the GB pbf and it worked fine. However, when I import the PBF into the db with the openstreetmap-carto style and use polytiles.py again, the tiles are created as before but without any large areas of water. e.g. the tile should look like:

but using carto and polytiles.py it looks like this:

I've run get-shapefiles.sh, carto is installed via 'npm install carto' and 'apt-get install node-carto' and I was wondering what steps I might have missed?


texts:
I have a Road Network shp file.(attributes include Names and Volumes )
Each road in my network consists of many links with the same name and different volumes.
I want to display Road Names and volumes as labels using sql query or python function to group by  names and have volume Average.

**
In normal SQL query I would type:

SELECT ROAD_NM, AVG(VOL) as VOL
FROM TABLE
GROUP BY ROAD_NM;


texts:
I wonder if you can help me resolve a problem.  I have given access to a WMS served up through Geoserver to various members of staff and it works fine on Arc 10 running on Windows XP. However I get an ArcMap drawing error when the WMS is called through Arc 10 running in a Windows 7 environment.

I get an ArcMAp Drawing Errors box appear with the following text within it:

One or more layers failed to draw:
Geoserver Web Map Service: WMS service exceptions:Service Error Description: Rendering request would use 1328087KB, whilst the maximum memory allowed is 1000000KB

The dilemma I have is I'm not sure if this is a memory allocation problem with Windows 7, Arc or Geoserver seeing as it runs without problem on XP.

If it is Geoserver, I'm presuming I would need to increase the JAI Maximum Memory; how do I do that?  The web isn't being very helpful I'm afraid.


texts:
I have a shapefile of survey control points (WGS 84) from the public domain provided by my local municipality.  I'm using a GeoXH with a zephyr antenna.  I have my coordinate system  set to WGS 1984.  I'm getting 4.5" accuracy (post processed) when I'm recording my test point.

I then take this data back to Trimble Pathfinder Office, post process it using a station about 3 miles from where the points were taken.  Pathfinder Office is configured to have WGS 1984 coordinates.  I then export this data to a shapefile (without a projection).

Using QGIS (2.18) I start a new project, select WGS 84 (EPSG: 4326), and load my city provided data.  I then load my test points (3 taken on tripod over one of the city's monuments).  I then load my shapefile, selecting WGS 84.  And my data is off by 5 ft (west, and a little north).

Any idea what I'm doing wrong?


texts: How to rotate rectangle which fits in a circle in leaflet

texts:
I have a circle (in code its called "cyc_166_1")  located at a certain angle with regard to another point.
I would like to draw a rectangle which both fits in that second circle (bounds within it), and has the same rotation angle.
How is it possible to do that?

var sainteAgathe_166_cyclic = L.circle([44.7105555555556, 6.84916666666667], 750, {fillColor:'none', color: '#404040', weight: 0.8, opacity: 1,}).addTo(map);
var length = 0.007
var angle = 55
var end_x = sainteAgathe_166_cyclic.getLatLng().lng + length * Math.cos(angle * Math.PI / 180)
var end_y = sainteAgathe_166_cyclic.getLatLng().lat + length * Math.sin(angle * Math.PI / 180)
var cyc_166_1 = L.circle([end_y, end_x], 90, {fillColor:'#800026', fillOpacity: 1, color: '#800026', weight: 0.5, opacity: 1}).addTo(map);


texts:
I am going to be using an AddJoin_management command in a python script at the end of several loops. To test, I used the AddJoin tool and use the below command in the Python window and made the parameters identical to what the code is(the parameters are correct for what I want to accomplish). When I run the AddJoin it crunches numbers for a while and then the whole program crashes and asks if I would like to send an error report.

arcpy.AddJoin_management("SignsToFinish", "SignID", "NewSigns", "sid", "KEEP_ALL")

That is the statement that is being used in both the script and in the AddJoin tool.

Why does this cause ArcMap to crash? Any ideas?


texts: Flat image output in TileMill SVG export

texts:
I have created a map in TileMill using OSM data (~500 Mb) and modifying the CartoCSS to fit my use case. I need to export it as an SVG with all the different objects being vectors etc.

However, when I use the built-in Export -> SVG function, the result is only a flat image of the area I defined in the viewfinder.

Except for some Font warnings during the export process, there are no errors.
How can I obtain a vector image of my map?

Here is my export window in TileMill, the resulting image opened in Inkscape (showing that it is indeed only a flat image and not a collection of vectors) and the terminal output (showing no error, only warning).


texts: Convert polygon data into extreme lon/lat coordinates using QGIS

texts:
I'm very new to QGIS. I have a shape file of census blocks and want to use the data to plot points within each boundary using Google Maps api.

I figure the best way to do this is to find the extreme north, south, east and west lat/long for each polygon (so later I can plot points within those boundaries).

I'm wondering:

Is there a way to convert the poloygons into lat/lon coordinates and
then find the 4 extreme boundary points?
Is there a better way to go about plotting points in boundaries that
I'm missing?


texts: Security alert dialogbox when launching arctoolbox

texts:
Every time I launch arctoolbox from arccatalog, the following dialog box appears:

Then, whatever I answer, it takes ages (~36s) before the box disappear and I can finally use arctoolbox. What could be done to remove that?


texts: I have a river raster that is composed of a range of canopy cover percentages by pixel. I want to rank which water basin's have the river with the best canopy cover. How do I do this?


texts: MapServer rendering at the lowest zoom levels

texts:
I am developing a web mapping application based on OpenLayers and using a WMS layer from MapServer (the source data for this layer is stored in PostGIS database) and I am facing the following problem: MapServer eliminates a lot of features at the lowest zoom levels.

Zoom Level 10:

Zoom Level 7:

Zoom Level 5:

Zoom Level 4:

Zoom Level 3:

I guess it is because MapServer doesn't display features with area less than one pixel. Is it possible to render such features at lowest zoom levels (for example as one pixel features)?

UPD.: I'm talking about orange features, green ones - it is part of OSM tiles.

Here is my mapfile:

MAP
SIZE 800 600
MAXSIZE 4096

IMAGECOLOR 255 255 255
IMAGETYPE PNG

OUTPUTFORMAT
    NAME "png"
    EXTENSION "png"
    MIMETYPE "image/png"
    DRIVER AGG/PNG
    IMAGEMODE RGBA
    FORMATOPTION "INTERLACE=OFF"
END

EXTENT -180 -90 180 90
PROJECTION
    "init=epsg:4326"
END

WEB
    METADATA
        wms_onlineresource "${request.route_url('mapserver.ows')}"
        wfs_onlineresource "${request.route_url('mapserver.ows')}"
        ows_title "mcx_pk"
        wms_enable_request "*"
        # ows_service_onlineresource  "${request.route_url('mapserver.ows')}"
        wms_srs "epsg:4326 epsg:3857 epsg:900913"
    END
END

LAYER
    NAME "fields"
    CONNECTION "user=pk password=******* dbname=pk host=***.info"
    CONNECTIONTYPE postgis
    PROCESSING "CLOSE_CONNECTION=DEFER"
    DATA "wkb_geometry from mask_f using unique ogc_fid using srid=4326"
    TYPE line
    DUMP TRUE
    TEMPLATE dummy.html
    PROJECTION
        "init=epsg:4326"
    END
    EXTENT -180 -90 180 90
    METADATA
        gml_include_items   "all"
        wms_include_items   "all"
    END
    CLASS
        STYLE
            COLOR 255 136 0
        END
    END
END


texts: Distribute ArcGIS services across multiple servers

texts:
What do you do when you have so many published services in ArcGIS that running even one instance of each service maxes out RAM/CPU?

We've grown to the point where we have more published services than our current server can run at one time, and expect to continue to grow considerably over the next few years. Even at the minimum one instance per service, we're running out of RAM.

Of course, doubling the RAM would get us by for probably another few months, but we will eventually need multiple servers, and I need to know how to work this out.

My understanding is that when you add a new server to an existing site, all of your published services will have at least one instance per service on all servers. Is there a way of distributing these services across multiple servers without having multiple 'sites'?

I konw that I could publish X amount of services to one server, and another batch to a separate server, but I'm assuming there's a better way to do this, and I'm hoping someone can shed some light on this.

Thanks!!!


texts: Licensing Issues while writing Community Project using ArcPy?

texts:
I'm creating a project for the Esri community, that uses some Basic ArcMap license, to create tools that require an Advanced (ArcInfo) license.

Is there a licensing issue to do something like this?

It's basically a Python Toolbox, that has some tools.


texts: Setting Output Parameter in ArcPy for Server Tool

texts:
I have published the following script tool (which is wrapped in a model, i.e. just the script tool and all of its parameters) to ArcGIS Server 10. The GP service succeeds when run, but I cannot determine how to access the final output of the GP service. The result is a RecordSet object, and I can view it in the arcgisjobs/job-id/ directory when viewing it on the server.

However, I cannot figure out how to access the result via the REST API. As you can see, the script tool doesn't have an output parameter. When viewing the GP service via the REST Directory there are only the two inputs. I have tried using SetParameterAsText and a Result object and adding an output parameter to both the ST and the model but without any success. So, the question is how do I modify the script and/or script tool/model properties to designate an output parameter so that it can be accessed via the REST API and subsequently consumed in a web application?

I have read through this related question but it has not helped.

The REST Endpoint: http://atlas.utc.edu/ArcGIS/rest/services/RRI/DE/GPServer/Decision%20Engine

The Python script:

#Import ArcPy site package module
import arcpy
from arcpy import env
from arcpy import Extent
from arcpy import Raster
from arcpy import RecordSet
from arcpy.sa import ZonalStatisticsAsTable
from time import gmtime, strftime

# Check out ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
# Set Environment Workspace
ws = env.workspace = WORKSPACE PATH

# Set Environment Raster Cell Size
cellSize = env.cellSize = 10
# Set Environment Processing Extent
processingExtent = Extent(602560.330672585, 3827530.65305605, 746490.330672585, 3966810.65305605)

environProcessingExtent = env.extent = processingExtent

# Strings representing rasters in the workspace
rasters = arcpy.GetParameterAsText(0).split(";")
# The dataset that defines the zones
inZoneData = arcpy.GetParameterAsText(1)

# Determine Zone Dataset
if inZoneData == "Block Group":
    zoneField = "Value"
    zoneData = ws + "\\ToolData\\Data.gdb\\BG_Ras"
elif inZoneData == "Census Tract":
    zoneField = "TRACTCE10"
    zoneData = "\\ToolData\\Data.gdb\\CT_Ras"
else:
    zoneField = ""

addedRaster = Raster(ws + "\\ToolData\\Data.gdb\\constantRaster")

# Process each raster string
for raster in rasters:
    addedRaster = addedRaster + Raster(ws + "\\ToolData\\Data.gdb\\" + raster)

# Execute Zonal Statistics
zonalTable = RecordSet(ZonalStatisticsAsTable(zoneData, zoneField, addedRaster, "%SCRATCHWORKSPACE%\\zontalTable", "DATA", "MAXIMUM")).save("%SCRATCHWORKSPACE%\\rszt")

Some relevant images?

Script Tool Parameters

Model Tool Parameters

Thanks.


texts:
I have a few NumPy arrays: lat, lon, and 6 data bands. I'd like to convert these to GeoTiff and then apply the gdalwarp function, based on a specific projection, to each band.

Is there a standard way to do this in Python?


texts: Can Geodatabase files from ArcMap be imported and used in Qgis?

texts:
I'm looking into Qgis and am doing some tests, but one thing I need to be able to do is use geodatabase files from ArcMap and ArcPad. Is that something that can be done and with relative ease?

I have been playing around with it for a few days and I haven't been able to figure it out.  I've also looked at several threads, but I haven't seen anything specific this topic. Any help or advice would be greatly appreciated.


texts:
I have two services, one a compact cached service the other an exploded cached service. Both are from the same data source, have the same publishing parameters (minus compact/exploded), and were published via ArcGIS Server 10.3.

In ArcMap, I make a WMS connection, add the exploded cached WMS service. It works fine.

Problem: I make a WMS connection, add the compact cached WMS service. The layer fails to draw, and eventually times out. In my ArcGIS Server manager log I get the following error:

I have tested a few things, and this is what I know:

The data frame default setting in ArcMap for WMS is EPSG: 4326, when you bring the service in
For the compact cached service, if I change the data frame to EPSG: 3857 the service draws
For the compact cached service, the service does NOT draw at EPSG: 4326, the default data frame setting
For the exploded cached service, the service draws if the data frame is in either EPSG: 4326 or 3857
I have tried to add EPSG: 4326 in the ArcGIS Server Administrator Directory > Edit Map Server > Service Properties, this did NOT work
The 'unspecified error' is not enough information to help
I do not get his issue if published via ArcGIS Server 10.2.2 (but I need to use 10.3)


texts: Software for building non-overlapping sets of polygons?

texts:
We have multiple populations of patients spread out over geographic distances.  For each population, we want users to easily define non-overlapping polygons.  The outer borders are likely to be U.S County boundaries, if that makes any difference.

Essentially, I'd like for users to visually build groups of non-overlapping polygons and have the application provide us with sets of coordinates for each polygon.  I would then be using the data in Google Maps.  We currently have a method to identify patient's location by lat/long, using Google Maps API if the user needs to approximate the location.

Web-based software is preferred, but I'd consider a desktop package if the interface was easy enough for our users with minimal training.


texts:
I'm trying to dissolve polygons based on a code number in the attribute table.

It all works well, with an exception of two polygons that fail to dissolve - they literally drop out from the output. It affects always the same two polygons. The dissolve process is successful, the layer is generated, just with two holes in it. ArcGIS does not give me any error message or anything.

I have tried messing with the source dbf file to see whether there are perhaps some funky characters in the code field, but it seems to be fine. The code is really just a simple number (12 and 707 specifically).

The only thing that seems to have perhaps some relevance is that both of those missing polygons happen to be parts of large cities, therefore they consist of higher number of small polygons to be dissolved (but still nothing that ArcGIS could not handle). Other than that I have really zero clue what could be wrong...

I'm using ArcGIS Desktop 9.3.


texts: Remove vertex from geometry (polygon) in OpenLayers 5?

texts:
I could able to remove the recent modified vertex of a feature using removePoint method of modify interaction, as well as using deleteCondition event of modify interaction. However, that is not what I am trying to achieve. The usecase for this this problem is as below.

I have modify interaction added on map. Now when user tries to modify a feature over map, I am highlighting the vertices, like this

so that user can directly go to these vertices and then do right click which will result in removal of vertex.

For this, I had to implement the contextmenu to handle the right click over map. But, now the challenge is to remove the received coordinate from the modifying geometry.

I tried to implement it as follows, with simple .splice() method on coordinate array of target geometry

 if(this.modifyInteraction.features_.getLength() && this.modifyInteraction.features_.getArray()[0]){
           //Feature being modified
           let modifyft = this.modifyInteraction.features_.getArray()[0];
           let coord = modifyft.getGeometry().getCoordinates()[0];
           let foundIndex;
           for(var k = 0; k < coord.length; k++){
              if(coord[k][0] == rightClickedVertex.getGeometry().getCoordinates()[0] &&
              coord[k][1] == rightClickedVertex.getGeometry().getCoordinates()[1]){
                 foundIndex = k;
                 break;
              }
           }
           if(foundIndex && foundIndex > -1){
              modifyft.getGeometry().getCoordinates()[0].splice(foundIndex, 1);
              //nothing changes
           }
        }

is there any other way to remove the vertex of geometry being edited?


texts:
I would like to implement the same flash feature command that ESRI uses when accessing the context menu on a table row.

I have used the following successfully to flash a selected feature:

IActiveView activeView = ArcMap.Document.ActiveView;
IFeatureIdentifyObj featIdentify = new FeatureIdentifyObject();
featIdentify.Feature = feature;
IIdentifyObj identify = featIdentify as IIdentifyObj;
identify.Flash(activeView.ScreenDisplay);

This code works, but does not have the 'crosshairs' additional effect that the ESRI flash command does.

I was looking possibly executing a command with ICommand.Execute, but could not find an ID for the Flash Command from the table context menu.

Anyone have success with implementing the same Flash Command that ESRI uses within custom code?    Does anyone know where to find the ID used for the Flash Command?


texts: Interacting with the Map after selecting a button/tool

texts:
I am using the ESRI Javascript API 2.0 and Dojo dijits.

Does anyone know how to wire up a dojo button (or toggle button) in a toolbar to a mouse click event on the map?

I want the user to have to select a tool/button before clicking on the map and NOT just have the map listen to any old mouse click event as this is too open i.e. I want to constrain the functionality to the selecting a tool from the UI.


texts: Would it be possible to find lat-long pairs for the center point of every section or even quarter section in a state? The accuracy of the point as the center is irrelevant, I would just need any point within the section or quarter section.


texts: Adding timeline play button in CARTO Builder?

texts:
I am using the new CARTO builder, trying to add a timeline with the time series widget. (something like this one: http://bl.ocks.org/csobier/cebdd47242d7ca98ec5e)

However, a bar chart shows up, which only allows for selection rather than playing the data overtime. I do see a play button in the thumbnail for type selection though, on the upper left corner.

I saw a lot of folks use Torque.js for building timelines, but what I am asking here is, is it possible to add a timeline play button inside the CARTO Builder?


texts: Feature class' default domain for a field is still locked to an older domain

texts:
I have a field within a feature class that I had assigned a coded text domain to, but then created a different coded text domain with updated information that I then pointed to as my new domain instead.

The new coded domain works fine for editing data and shows no sign of the original domain, but when I go to delete my old domain I no longer use from the file geodatabase it tells me 'domain is used as default domain'. I get this error whether I go through the domain options within the feature class in ArcCatalog or through the delete domain tool in ArcToolbox/Data Management/Domains.

ESRI File geodatabase was created and is being edited in Arc 10.0. on the same machine by the same user.

How would I go about removing this domain error or do I need to duplicate the FGDB and hope it removed any corrupted files?


texts:
I have a layer with a lot of polygons. Now I want to filter all polygons that are very narrow.
For example, the polygons show living space. But I'm only interested in large living spaces. Function $area does not help here, as there are also very long and narrow surfaces.
Does anyone know a function in QGIS that comes to this output via filtering x and y?


texts: Generate ArcGIS spatial bookmarks from feature envelope

texts: Curious if anyone has some code that would rip through a feature class and generate a spatial bookmark for each feature's envelope. There can be a field in the feature class to hold the bookmark's name.


texts: How are OSM tags mapped to fields in Mapbox Streets vector tiles?

texts:
I recently started experimenting with Mapbox Studio and have encountered some discrepancies between feature tags in OSM and the corresponding feature fields in the vector tiles from Mapbox Streets.  The Mapbox Streets documentation for the roads layer states that the type field for most roads will be the value of the OSM highway tag.  However, it appears that in OSM the following feature has the value 'path' for the highway tag while the attached image showing a query of the feature with the same id in the Mapbox Streets vector tiles shows it type value as 'piste' rather than 'path'.  I encountered several instances of this that were giving me some unexpected styling results in one region.

How are OSM tags mapped to vector tile fields?


texts: Problems when use GeoPandas for Itersection between two Lines?

texts:
How to a get the points of intersection between two geodataframe with different Z?
In dataframe I Have just the name of line(index) and Geometry LineString Z.

I try this code but, always the intersections give the same Z value

I have two geodataframe df_cross and df_lines

def load_xyz(file_path, file_name, spatial_reference, sep=';', header=0):
df = pd.read_csv(file_path+'\\'+file_name, header=header, sep=sep,
                 names=['Line','X','Y','Z'])

spatial_reference = {'init': 'epsg:4326'}

gdf = gpd.GeoDataFrame(df.values,
                       crs=spatial_reference,
                       geometry=[Point(xyz) for xyz in zip(df['X'], df['Y'], df['Z'])],
                       columns=['Line','X','Y','Z'])

gdf = gdf.groupby(['Line'])['geometry'].apply(lambda x: LineString(x.tolist()))

return gpd.GeoDataFrame(gdf)

spatial_reference = {'init': 'epsg:32724'}

file_path = r'C:\Users\VN\Desktop'
file_name = r'xyz.csv'
df_lines = load_xyz(file_path, file_name, spatial_reference)

file_path = r'C:\Users\VN\Desktop'
file_name = r'xyz_cross.csv'
df_cross = load_xyz(file_path, file_name, spatial_reference)

for i in range(len(df_lines)):
    for j in range(len(df_cross)):
        if(df_lines.iloc[i]['geometry'].intersects(df_cross.iloc[j]['geometry'])):

            df_lines.iloc[i].geometry.intersection(df_cross.iloc[j].geometry).z

            df_cross.iloc[j].geometry.intersection(df_lines.iloc[i].geometry).z

Does anyone know the way to solve this problem?


texts: How to indicate a change in transport mode for directions in a multimodal dataset

texts:
When I solve a route and generate directions for my multimodal dataset, the end product is not very user friendly:

The pedestrian started walking along Main Street and then boarded a bus to travel along the various segments of the bus route. Is there a way to change the text so that it makes more sense, so something like:

Walk northwest along Main St
Board the bus at Stop 1
...etc

The connectivity is set up so that the pedestrian can only get on the bus at official bus stops.


texts:
I'm working on a series of maps and whenever I select a symbol, or pan, or anything my whole map has to redraw.  It is a very slow redraw.  I generally have to wait 5-10 seconds, even if I am just selecting one symbol to move.  Also, it doesn't draw everything, so I think I have moved all the grills from behind trees, but inevitably there is always one that does not draw.  The imagery I am using as a basemap is a very large raster, though even when I turn it off, it still takes a good chunk of time redrawing.

Does anyone know what the problem might be?  I thought it might be the graphics card but I'm not sure how to tell or if it might be another issue.

Here are my computer stats:
Operating System: Windows10 Enterprise 64-bit
Processor: Intel(R) Core(TM) i5-4590 CPU @ 3.30GHz (4 CPUs), ~3.3GHz
Memory: 8000MB RAM
For both monitors:
   Name: Interl(R) HD Graphics 4600
   Device Type: Full Display Device
   Approx. Total Memory: 4167 MB
   Current Display Mode: 1920 x 1080 (32 bit) (60 Hz)


texts:
I am still new to this, so I apologise if I don't make a lot of sense!

I am using QGIS 2.4.0 Chugiak.

I am working on a project set to the GDA94/MGA Zone 50 coordinate system with 'on the fly transformation' enabled. I have added a shapefile that was created in ArcGis and set the coordinate system to the same as the projects (GDA94/MGA Zone 50). However, the shapefile is projecting at a much larger scale and A LOT further away than it should be. I have tried changing the coordinate system of the shapefile and reprojecting the shapefile, but haven't had any luck.

Any chance someone may be able to help?


texts: Issue with initializing ArcObjects Interface and classes

texts:
I was going through ArcObjects example on following linkWalkthrough_Building_custom_UI_elements .I followed same steps but I was getting error while initializing objects for following code.

ESRI.ArcGIS.Display.IRubberBand rubberBand = new ESRI.ArcGIS.Display.RubberLineClass()
ISimpleMarkerSymbol simpleMarkerSymbol = new ESRI.ArcGIS.Display.SimpleMarkerSymbolClass()
IMarkerElement markerElement = new MarkerElementClass()

etc I was not able to find exact reason why I was getting but may be due to version of ArcObjects.I am using .Net 4.0,Visual studio 2010, ArcObjects 10.2.
I got following type of error

The type 'ESRI.ArcGIS.Display.RubberLineClass' has no constructors defined.
Interop type 'ESRI.ArcGIS.Display.SimpleMarkerSymbolClass' cannot be embedded. Use the applicable interface instead.


texts: Visualizing all Sentinel-1 Images in collection using Google Earth Engine

texts:
Here is the code. There are 81 sentinel images but I am only seeing the first one (or the last one maybe), how can I visualize other images. Let's say I want to see the 2nd one ?
// Importing S-1 Images
var S1 = ee.ImageCollection("COPERNICUS/S1_GRD");

// Importing Geometry (area of Interest)
geometry =
ee.Geometry.Polygon(
    [[[74.38065032641738, 12.778968253582836],
      [74.38065032641738, 7.612293582692355],
      [77.74246673266738, 7.612293582692355],
      [77.74246673266738, 12.778968253582836]]], null, false);

// Define Time Interval
var start_date = ee.Date('2018-07-21');
var end_date = start_date.advance(3, 'month');
var date_filter = ee.Filter.date(start_date, end_date);

//Display of Image
var collection2 = S1
    .filter(date_filter)
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filterBounds(geometry);
var vis_params2 = {bands:'VV', min:-25, max:5};
var label2 = 'SAR Backscatter';

// Visualizing the map
Map.addLayer(collection2, vis_params2, 'secondimg')
Map.centerObject(geometry);


texts: Integrating use of surface distance when calculating cumulative cost-surface based on walking pace?

texts:
Disclaimer

This question is somewhat related to a number of earlier questions of mine:

Accumulated cost surface (Tobler's hiking function) ignoring slope using R gdistance?;
Calculating accumulated walking pace from starting location outwards using R?;
Working out travel time from gdistance's conductance values in R?,

each representing a step in my struggle to implement in R the calculation of isochrones around a given location, using the Tobler's hiking function. In what follows, I am going to ask something new, albeit it will entail making reference to what I have achieved in an attempt to address the problem I described in my earlier threads. Long story short: I hope this new question will not be flagged as possible duplicate.

Background

As said above, I was trying to calculate isochrones around a start location using the Tobler's hiking function (solved for time). I believe I eventually got something right, building on the gdistance package's documentation (https://mran.microsoft.com/snapshot/2014-11-17/web/packages/gdistance/vignettes/gdistance.pdf).

Below is a reproducible example, before I ask my new question:

First, load the gdistance library and a sample dataset (DTM):

library(gdistance)
dtm <- raster(system.file("external/maungawhau.grd", package="gdistance"))

Secondly, calculate the terrain slope:

altDiff <- function(x){x[2] - x[1]}
hd <- transition(dtm, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd)

Thirdly, apply the Tobler's function to only those slope cells that are adjacents:

adj <- adjacent(dtm, cells=1:ncell(dtm), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- (6 * exp(3.5 * abs(slope[adj] + 0.05)))^-1 #s/m

In the above step, I used the reciprocal of the speed in m/s. Rationale: if we use speed, the final accumulated values will be 1/travel time (according to gdistance pdf guide; link above); therefore, I use the reciprocal of speed to eventually get travel time/1; travel time will be in seconds, and will require to be divided by 3600 to get hours.

Fourthly, geocorrection of the speed values (well, we should call "pace" not "speed"):

Conductance <- geoCorrection(speed)

Finally, accumulating the speed (i.e., pace) from a location:

A <- c(2667670, 6479000)
accum_sec <- accCost(Conductance, A)
accum_hour <- accum_sec/3600

What we get using rasterVis is the image further below:

levels <- seq(min(values(accum_hour)), max(values(accum_hour)), 0.15)
rasterVis::levelplot(accum_hour, par.settings = RdBuTheme, margin=FALSE, contour=TRUE, at = levels, lty=3, labels = list(cex = 0.7), label.style = 'align')

Given the background, here is the question

While the result looks quite good, in order to best mimic part of what is done by arcGIS's Path Distance tool (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-the-path-distance-tools-work.htm), I am wondering if is it possible to take into account surface distance in the calculation of the accumulated cost-surface, not the planar distance?

In other words, when passing from cell A to B, the cost should be (as customary) the average of A and B (since the traveller experiences half of the cost of the origin cell and half of the cost of the destination cell), but multiplied by the surface (not planar) distance (which could be quite easily derived using the Pytagorean theorem knowing the elevation difference between A and B, and their distance).

While in theory the above sounds easy, I cannot figure out how to actually integrate that to the workflow I have described earlier in this thread.


texts: How do you clip a gnomonic or stereographic azimuthal projection?

texts:
I have polygons comprised of lat/lon coordinates and I want to project them onto a plane that's appropriate for polygon operations (such as checking for intersections). In some cases, the polygons that I want to project might fall outside the valid area that can be projected (such as gnomonic being confined to less than a hemisphere).

Given an arbitrary center lon and lat, how can I clip polygons that are partially outside the valid range so I can project them? In the attached picture, note that the polygons and graticules are clipped.

The way I do a gnomonic projection currently is directly project coordinates in cartesian (ECEF) out from the center of the Earth onto a tangent plane. But I obviously can't project any points or lines where ray to a coordinate is greater than or equal to 90 degrees from the center lon and lat.

I'm looking for a method description as opposed to a pre-existing software package.


texts: Remove all joins from a layer using pyqgis

texts:
The title says it all. I want to remove all joins from a layer.

After examining the QGIS API documentation and playing a bit with the QgsVectorLayer class I've figured out how to do it but my doubt is whether there's a better way to do this same thing. I haven't found a single example using QgsVectorLayer.vectorJoins nor anything related to identifying which tables are joined to a layer so I don't know if this is actually the way to go.

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    if layer.name() == "Fire Halls":
        fh_lyr = layer

joinsInfo = fh_lyr.vectorJoins()
for joinInfo in joinsInfo:
    ok = fh_lyr.removeJoin(joinInfo.joinLayerId)
    if ok:
        progress.setInfo(u'Join removed successfully')
    else:
        msg = "Couldn't remove the layer join"
        raise GeoAlgorithmExecutionException(msg)


texts:
I am using ArcGIS Java Script API version 2.8 & using query layer in ArcGIS server.

The data is coming from sql server and most of the data is in Arabic. I am sorting the data according to  the particular map is displaying.

I have created the query task and the result will be displayed in a combo-box but when I select Arabic record, the map is displaying blank. If I select English record then its displaying correctly. For details please check the screenshot.

The column type is nvarchar in sql server and in query task, I am using the where condition as follows but still getting an error:

query.where = "STATE <> N''"  (In sql server query works if I use N'<Arabic text>')

What I tried:

Changed the query parameter.
It works perfectly for English names.
Added UTF 8 encoding.
Upgraded API to 3.2 (latest one)
In Sql server where query works successfully.

Still getting an error. Any help will be great!!! thanks in advance : )

@All

I confirmed through REST API and through API that I was able to query to the database. Please check the new screenshots

So as far as the logic, the information is coming but its not displaying on the map (its interesting), also I have checked parameters in the map service and still researching.


texts:
The new Mapbox Studio seems to have a button for downloading a style:

but not a dataset:

So, is there actually a way to export a whole dataset?

This may seem like a strange thing to do, but since the editor is now cloud-based, you don't necessarily have the source files close at hand.

The only way to download the data again seems to be to request a very large number of vector tiles, and somehow stitch them back together.


texts:
I would like to use the distance matrix tool (or a similar function) to look at distances between sites (displayed as point data) and large rivers (displayed as linear).

I have been following this nearest neighbour tutorial until now
http://www.qgistutorials.com/en/docs/nearest_neighbor_analysis.html

But the distance matrix tool only allows for point data. Apart from the fact my input layer (European rivers) is linear, the initial aim of this research is pretty basic. I want to find the mean distance between my sites and major rivers (later i want to examine the results in terms of travel times using slope data but that's a diffrent question I guess).

Is there a precedent for this kind of research or am I on the wrong track entirely? Am using QGIS 2.16.


texts:
I am trying to extract values from a landuse raster layer to a point shapefile. However, when I use either Extract Multi Values to Points or the Extract Values to Points tool, it returns pixel values that do not resemble the land use classification. For example, the land use code system looks like this:

1.1 Nature conservation
1.1.1 Strict nature reserves
1.1.2 Wilderness area
etc.

How can I extract this information from the raster into a point feature class?


texts: In QGIS 1.8, I cannot find the way to assign a field to 'map tips'. There seems to be no way in 'properties' of the layer to do so.


texts: How to setup a new script, and parameters in ArcMap?

texts:
I found a python script (Drawing parallel lines inside polygons (Well Paths)?) I would like to use in an analysis. However, I don't understand how to add the script posted by @phloem to ArcMap.

Following an an older tutorial, for ArcMap 9.2 (http://webhelp.esri.com/arcgisdesktop/9.2/tutorials/GP_3_3.htm) assisted me in the creation of a script within a toolbox. But, after running the new script, the dialogue box says "no parameters".

How do I setup the new script, along with the relevant parameters necessary to execute the script as a tool? I am using ArcMap 10.3.1


texts: How to add ITable (Table object) to current map (or document)?

texts:
I created an in-memory results ITable and I need to add it to the current ArcMap document/map to see the generated data.

How can I do it?


texts: How many satellites are actually used to get position?

texts:
NMEA 0183 Version 4.10 defines that there cannot be more than 12 satellites in GGA message:

Number of satellites used (range: 0-12)

But there can be more active satellites defined in GSA:

satellites used for navigation, and DOP values

So how many satellites are actually used to calculate position?
See example:

$GNGGA,155000,4700,N,03201.46451,E,1,12,0.78,9.5,M,27.7,M,,*4E
$GNGSA,A,3,21,18,20,16,10,15,27,08,32,,,,1.27,0.78,1.01*11
$GNGSA,A,3, 68,83,76,69,85,,,,,,,,1.27,0.78,1.01*11

12 satellites in GGA
14 (9+5) satellites in GSA


texts: Displaying both data in different fields in QGIS

texts:
I have a Shapefile that is automatically loaded in QGIS. While loading the data the coordinates x & y are automatically displayed.
In my attribute table I have lat_1, long_1 and lat_2, long_2 in the same row. How can I display both at the same time? Is there any solution for that in QGIS?
lat_1, long_1 is the main location. lat_2, long_2 is entry point location.
That means in attribute table field names should be unique right. So it will display only x and y. But in my attribute table there are entry points (lat, long) they are also coordinates. So I want the both x & y (main and entry) displayed. Added screenshot is the attribute table format


texts:
As a follow up to Correctly setting up CRS systems for the NUTS (Eurostat) geographies and Doogal's data? I have spatially matched a set of files, one with lat/long coordinates for postcodes sourced from Doogal's website. The over function available through the sp in r package seems to generating plausible results. After assigining crs somehow arbitrarily I'm interested in verifying that the obtained spatial match is accurate.

Are there any methods I could employ to ensure that the obtained spatial match is error-free?

Reproducible example

# Download ----------------------------------------------------------------

if (!dir.exists("NUTS_2013_01M_SH")) {
    tmpzip <- tempfile(fileext = ".zip", pattern = "estat")
    download.file(url = "http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2013_01M_SH.zip",
                  destfile = tmpzip)
}
if (!dir.exists("NUTS_2013_01M_SH")) {
    unzip(zipfile = tmpzip,
          overwrite = TRUE,
          exdir = getwd())

}

# Import ------------------------------------------------------------------

shpsEurs <- maptools::readShapePoly(fn = file.path(getwd(),
                                                   "NUTS_2013_01M_SH",
                                                   "data",
                                                   "NUTS_RG_01M_2013"),
                                    verbose = TRUE)

# Projections -------------------------------------------------------------

# Create projections
wgs84_crs <- sp::CRS("+init=epsg:4326")
## Assign
sp::proj4string(shpsEurs) <- wgs84_crs

# Source PCs --------------------------------------------------------------

if (!file.exists("Scotland postcodes.csv")) {
    download.file(url = "https://www.doogal.co.uk/UKPostcodesCSV.ashx?country=Scotland",
                  destfile = "Scotland postcodes.csv")
}

dtaPCS <- read.csv("Scotland postcodes.csv", stringsAsFactors = FALSE)

sp::coordinates(dtaPCS) <- ~Longitude+Latitude
sp::proj4string(dtaPCS) <- wgs84_crs

# Match -------------------------------------------------------------------

res1 <- sp::over(x = dtaPCS,
                 y = shpsEurs)

Results

>> head(res1)
  NUTS_ID STAT_LEVL_ SHAPE_AREA SHAPE_LEN
1     UKM          1   11.62568  127.5792
2     UKM          1   11.62568  127.5792
3     UKM          1   11.62568  127.5792
4     UKM          1   11.62568  127.5792
5     UKM          1   11.62568  127.5792
6     UKM          1   11.62568  127.5792


texts: Selecting random points based on the distance between them

texts:
I have to select a random set of points from a point shapefile based on the following rule:

"no two points can be within 1km of each other"

Here is an example:

Ideally, I'd like to do this through the algorithms already present in QGIS, but if I having to do that through scripts it should not be a problem.


texts: Taking Union of several geometries in GEOS Python?

texts:
I am using the GEOS/OGR library in Python.  I have files full of geometries that overlap, and I want to convert these to files of geometries that don't overlap.

The way I have been approaching this is to just read in all of the geometries from a file, and then looping over each geometry and taking the Union() of that with all of the previous geometries.  So something like:

def dedup(geometries):
   """Return a geometry that is the union of all geometries."""
   if not geometries:  return None
   current_union = geometries[0].Clone()
   for g in geometries:
      current_union = current_union.Union(g).Clone()
   return current_union

Most of the time this appears to work, but occasionally Union() will return None.  My debugging reveals that under these cases:

IsValid() returns true for both g and current_union
Intersects() return true between g and current_union
GetGeometryCount() for current_union typically returns a number over 100.
I suspect the error might be some kind of memory issue because it seems like the error "moves" depending on context.  Ie if I say run this on files 1-100, it might crash on file 20, but if I say run it on files 20-100, it might crash on file 40.  I never use Destroy() on any geometries. Task Manager does not show me running out of memory though.

Anyway, maybe there is an easier way to do this?  Or else does anyone have any idea why it is crashing here?


texts: How can I edit lines containing M values with the ArcGIS JS API?

texts:
I am working with the ArcGIS JavaScript API and want to move and edit the vertices of a polyline analog to the following sample: http://help.arcgis.com/en/webapi/javascript/arcgis/demos/ed/ed_toolbar.html

This works fine as long as the polyline does not contain M Values. If it does contain M Values I get the following error in firebug.
Errorcode:-2147221399, description: „Geometry does not have M Value(s), hence cannot update the feature."

What am I doing wrong?

Here’s my code:

function initEditing(layers) {
    var myMap = this;
    var firePerimeterFL = myMap.getLayer("firePerimeterFL");

    var editToolbar = new esri.toolbars.Edit(myMap);
    dojo.connect(editToolbar, "onDeactivate", function(tool,graphic,info) {
      if(info.isModified){
        firePerimeterFL.applyEdits(null, [graphic], null);
      }
    });

    var editingEnabled = false;
    dojo.connect(firePerimeterFL, "onDblClick", function(evt) {
      dojo.stopEvent(evt);
      if (editingEnabled === false) {
        editingEnabled = true;
        editToolbar.activate(esri.toolbars.Edit.EDIT_VERTICES , evt.graphic);
      } else {
        editToolbar.deactivate();
        editingEnabled = false;
      }
    });
  }


texts: Trouble shooting freezing/stalling issues in QGIS

texts:
I've been using QGIS for more than a year now (currently using 2.8.2 on MacOS 10.10.3, but also in Windows XP via VirtualBox for troubleshooting) and the application has been freezing and stalling quite a lot on me as I've been trying to do some analysis over the past several weeks such that I'm nearly at a standstill.

The sort of work I'm doing is relatively mundane doing analysis using a set of a few hundred data points (in CSV) and a series of shapefiles (Scottish Urban-Rural data, Deprivation, etc.). I've struggled to isolate the source of the failure and have tried a variety of things (converted CSV to shapefiles, messed around with using different CRS so that all are matching, have tidied up my CSV data, reimported data and layers into a fresh QGIS file, etc.), but can't seem to discover what is causing my issues.

I expect that having a better grasp on good practices for troubleshooting QGIS might help me to isolate what is going on, and so I'm wondering if a savvy user can recommend some good practices for troubleshooting this kind of issue in QGIS (when analysis operations are freezing or stalling). Is there a way to view progress in a console window? Or maybe there are some logfiles I should be checking for specific errors?


texts: GeoServer NetCDF plugin, read grid_mapping_name="vertical_perspective"

texts:
I am currently using NetCDF plugin for GeoServer. I would like to know if there is any solution to be able to read a NetCDF with the attribute grid_mapping_name="vertical_perspective".
While debugging I found out in org.geotools.coverage.io.netcdf.crs.NetCDFProjection, that only few projections are currently supported (code available here : https://github.com/geotools/geotools/blob/fe5eac60bcbaabc50af40b75de7bae6f320f09e9/modules/plugin/coverage-multidim/netcdf/src/main/java/org/geotools/coverage/io/netcdf/crs/NetCDFProjection.java#L236)
So I would like to know if any of you dealt with "vertical_perspective" projection while using Geoserver NetCDF plugin? And how did you manage to make it work? Maybe one solution is to override NetCDFProjection class and add "vertical_perspective" projection?


texts: Can I put a button that automatically inputs GPS coordinates into a google form? (or other spreadsheet populator)

texts: I am currently creating a Google Form and part of the questionnaire requires the user's current GPS location. Is it possible for a button to be able to do this and automatically populate the field? Or is there another way to do this without Google Forms but still can export to Excel?


texts: Using ArcGIS Field Calculator and Python Parser to add zeros and ones?

texts:
I have a lot of rows of data, and I want to use the field calculator to give a uniform assignment of

0
1
0
1
0
1
...

but I don't know how to implement it.

Can someone help me?


texts:
I have a shapefile of Ireland, and some xy coordinates plotted on top of it.

Is it possible to calculate distance to the coastline, for each XY?

I am an ArcGis beginner.


texts: QGIS r.watershed ERROR: Raster map or group not found

texts:
Yesterday I was able to run the r.watershed tool multiple times without any issues, today it's failing every time with "ERROR: Raster map or group not found".  I am running the same DEM with the same parameters.  This is in QGIS 2.18 but I get a similar error when attempting the same in QGIS 3.

I completely uninstalled and re-installed QGIS but this did not help.  Any solutions on how to fix this?

EDIT:
I decided to run this in GRASS itself and there are no issues.


texts: Web Mapping with Select by Rectangle and Modify Symbology?

texts:
I want to create a free/simple web map where users can click & drag to select polygons in the map, then once selected, choose a color for all selected polygons. The map should also maintain these changes even after the user is done (not just browser-side styling).

For example: A simple web map containing polygons of the 50 US states is created. Dave opens the page, clicks & drags to select Texas and New Mexico, and chooses the color red, which colors both states red at once. Later, John opens the page and sees that Texas & NM are red, and then clicks and drags to select New Mexico, Arizona and California and changes their color to blue at once.

I'm experienced in using various mapping libraries (ArcGIS JS API, Google Maps API, Leaflet, etc), so it doesn't really matter which I use, but it just needs to be as simple as possible and free (e.g. no ArcGIS Server). This isn't a commercial application - purely recreational, but I can't seem to think of a simple solution for this.


texts: How to reload a georeferenced image in QGIS 2:16?

texts:
I georeferenced an image in QGIS 2:14 and these are the resulting files:

RawBingHybr17.bmp
RawBingHybr17_georef.tif
RawBingHybr17.GCP.points
RawBingHybr17_GCPpoints.points

I´m note sure if they are correct because I'm not an experienced user and has done this for some time ago.

I delete the layer from the project. In the QGIS Browser I saw the layer SRC as EPSG:4326 (WGS 84)

Now I want to insert it in another project in QGIS 2.16, with SRC EPSG:3857 (WGS 84 / Pseudo Mercator). When I insert the image on a raster layer, I see the image, but is away from the correct coordinates.

I will need to start over and georefence it again?

EDITED LATER:

I loaded another tif file georeferenced and was properly in place.

SOLVED:

In Georeferencer I loaded this tif georeferenced file and the points file. When running the Georeferencer, this one was put in its place..


texts:
Has anyone had experience successfully deploying GeoDjango with MongoDB? Is it possible? If so, how is it as simple as as installing the non-relational fork of Django and setting MongoDB as the backend?

I like the scaling capabilities of MongoDB and the ability to make schema changes without having to migrate my apps, but have had difficulty finding documentation on setting GeoDjango with non-relational backends.

My geospatial queries largely consist of checking for intersections and containment of geometries.


texts: Verifying Many Hyperlinks in an Attribute table at Once

texts:
I have been using ArcMap for a few years, but I have zero experience with Python or ArcPy. I have been looking for a way to check the links in my attribute table to see if they are broken, and by that I mean, which ones open to a 404-file not found page in a web browser . All my hyperlinks in the attribute table look like https://GIS.com/AsBuilts/PDF-NAME-HERE.pdf. I have had some exposure to other languages so I've been familiarizing myself with the code from this post: Using ArcPy to verify paths of Hyperlinks?. I got it working, but it always outputs all of rows in my attribute table to the broken links text file.

Which would mean that it thinks that all of my hyperlinks are broken?

I know for sure that there are some broken links and some that are good. I have linked them to the hyperlink tool in Arc and they open just fine.

I am not sure where to go from here, any suggestions of things to check?

Here's my code:

import arcpy
import os

fc = r"Z:\workplace\AsbuiltPoly_copy.shp"
fields = ["ASBUILT_RE"]

with arcpy.da.SearchCursor (fc,fields) as cursor:
    for row in cursor:
        Roll = str(row[0])
        if os.path.lexists(Roll):
            pass
        else:
            logfile = open(r"Z:\workplace\BrokenLinks.txt","a")
            logtext = Roll + os.linesep
            logfile.write(logtext)
            logfile.close()

I cannot download a link checker onto this computer cause I am not the admin :/ so I cannot go that route. But if you know anyone that are good, send them my way.


texts: Setting global CSS and Javascript in .KML

texts: I've tried googling this question, but I was wondering if it is possible to declare common styles (CSS) and Javascript functions that are being referenced within the  tag of multiple placemarks?


texts: Could SQL performance be increased by storing rotations of the data?

texts:
I'm working on a project where this kind of SQL Psuedocode would apply:

SELECT * FROM objects WHERE (distance

My inputs are latitude and longitude, and my database holds latitude and longitude.

One solution that I have seen is to use a command like the following to produce a square shaped output

SELECT * FROM objects WHERE (latitude BETWEEN center-distance AND center+distance) AND (longitude BETWEEN center-distance AND center+distance)

Then have the program shave off the corners of the square.

Does anyone know if it would perform well to make 2 extra rotations of the data and search by that? The output would be a 12-sided object. It may look something like this

SELECT * FROM objects WHERE(latitude BETWEEN a1 AND a2) AND (longitude BETWEEN b1 AND b2) AND (latitude_r1 BETWEEN a1_r1 AND a2_r1) AND (longitude_r1 BETWEEN b1_r1 AND b2_r1) AND (latitude_r2 BETWEEN a1_r2 AND a2_r2) AND (longitude_r2 BETWEEN b1_r2 AND b2_r2)

So, for example, no added rotations would make a square, 1 extra rotation would make an octogon (8 sided), 2 extra rotations 12, 3 extra rotations 16, meaning each additional rotation would increase accuracy, but I don't know if the additional indexes would cost more or less than the gains per a rotation.

So, assuming that creating the rotated data set and creating the rotated search terms were free, would this be a viable method to estimate objects within X distance? In my situation, performance is more important than 99% accuracy.


texts: QGIS 3.18 - Save form field value into variable

texts:
I'm trying to find a way to store a value from a QGIS form field into a QGIS project variable. Is it something possible ?
I have a postgreSQL table and a QGIS form made from the fields of that table. What I'm trying to do is once the form has been submitted, the value of the id field is stored into a QGIS project variable to be used in another form.
If that's not possible using the default QGIS forms functionality, can it be implemented using python ? I didn't find any documentation about working with forms and python.


texts: Visualizing and organizing traffic flow data for an intersection

texts:
I'm looking for advice on how to organize data that describes traffic flows through several unconnected intersections over several years. I'd like to learn about the most basic representation of this data, that would lend itself to future more complex analysis. In this case the traffic flow is of bicycles, so the numbers are quite low.

I read through this question and all the answers. glennon provided several useful links in his answer, but I think it that such an approach is more sophisticated than required here.

I'm also interested in compelling ways of showing these flows over time and would appreciate any suggestions.


texts: How to find all cells within a neighborhood with slopes towards the center?

texts:
I would like to iteratively go through each cell in a raster in arcGIS that meets criteria (already established as per below description) and find the percentage of cells in a several km neighborhood with slopes pointing towards/away from the central cell (before moving onto the next cell, which would then become the central cell).

For a given raster of a region (let's say Europe so you have the idea of the extent of the dataset) I have created a map of all slopes above a certain degree.
From this raster, I am calculating a ruggedness variable that applies to those slopes (effects things like wind speed) using focal statistics: percentage of cells in the specified neighborhood of each cell that is above that slope. But, the ruggedness variable does not truly capture this when it includes slopes that do not effect the wind.

An example would be that when looking in the direction due north of a specified cell, within the neighborhood, a slope is oriented east-west. While this has some effect on the wind, it is not the specific effect I am looking for. I would want to only look at slopes oriented north, northwest, northeast (to a certain limit) and south, southeast/west that lie north (or south) of the central cell, and then also for cells east or west of the cell, only slopes oriented east/west/etc.

I have the idea of using Aspect on the slopes meeting my criteria, which would give the angle of each cell. I thought about using Euclidean direction and iterating through each cell, defining that cell as the source, and then applying an equation to cells in the neighborhood (a few km) to determine if the aspect and euclidean direction match. However, in order to do this I would need to basically have a custom focal statistics tool.

Has anyone done something similar or do you have any suggestions?


texts:
I have a layer in PostGIS in epsg:4326. Put in in MapServer:

MAP
IMAGETYPE PNG
EXTENT 60.0 56.0 61.2 57.4701
SIZE 800 600
OUTPUTFORMAT
    NAME "png"
    DRIVER        GD/PNG
    IMAGEMODE     PC256
    TRANSPARENT   ON
END
    PROJECTION
           "init=epsg:3857"
    END
LAYER
    NAME "boundary"
    CONNECTIONTYPE postgis
    CONNECTION "user=postgres password=postgres dbname=ekb_1 host=localhost"
    DATA "the_geom from bige"
    TYPE POLYGON
    CLASS
            COLOR 200 120 123
    END
END

END

Now i add it to leaflet map:

            var map = new L.Map('map',{
                center: new L.LatLng(56.8391677245,60.7011796115)
                ,minZoom: 12
                ,maxZoom: 20
                ,zoom: 14
                ,zoomAnimation: true,
                measureControl: true,
            });

            var nexrad = L.tileLayer.wms('http://127.0.0.1/map/?map=/srv/www/localhost/public_html/map/ekb_test.map&mode=map', {
           layers:'boundary',
            format: 'image/png',
            transparent: true,
            attribution: "Weather data © 2012 IEM Nexrad"
        });
        map.addLayer(nexrad);

And see multyple layer images on map. What can be wrong?

UPDATE

There is my  /usr/lib/cgi-bin/mapserv -v output:

MapServer version 6.2.1 OUTPUT=GIF OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ
 SUPPORTS=GD SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO
SUPPORTS=SVG_SYMBOLS SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER
SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT
SUPPORTS=WCS_SERVER SUPPORTS=SOS_SERVER SUPPORTS=FASTCGI SUPPORTS=THREADS
SUPPORTS=GEOS INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE


texts: How to Add Field to a Shapefile Using DotSpatial

texts:
I tried a lots on Google to find a sample of showing how to add a Field to attribute table of an existing shapefile. For example I have a Shapefile at

  C://data/Streets.shp

and need to add two field L_CITY and R_CITY both text and 50 characters limit. Can you please let me know how to do this in DotSpatial?


texts:
I have a series of bandmath processes that need to be run via envi and that I am doing over and over again. I have no experience with IDL but am trying to automate using it. I want to be able to choose input files via dialog box, or even just be able to pass the result of a bandmath process to the next process. I am able to do this automation using python and arc but don't seem to get the correct results, hence IDL and ENVI.

IDL documentation seems to indicate that that in IDL the bands need to be in the same raster whereas in ENVI I never need to do that. Is this the case? Am I missing a step to get my code to work?

Barring specifics, are there better tutorials out there than what's available in the code documentation?

My code

; Start ENVI
e = ENVI()

; Open a data file
File1 = DIALOG_PICKFILE(/READ)
Raster1 = e.OpenRaster(File1)

; Determine an output file
OutFile1 = e.GetTemporaryFilename()

File2 = DIALOG_PICKFILE(/READ)
Raster2 = e.OpenRaster(File2)

Outfile2 = e.GetTemporaryFilename()
; Return a file ID
fid = ENVIRasterToFID(Raster1)

Outfile3 = e.GetTemporaryFilename()

; Add the first two bands
ENVI_File_Query, fid, DIMS=dims
ENVI_Doit, 'Math_Doit', $
  FID = [fid, fid], $
  DIMS = dims, $
  POS = [0,1], $
  EXP = 'b1 + b2', $
  OUT_NAME = Outfile3

END


texts: Allowing all hosts on map in Mapfish Print 2.1.2?

texts:
My web mapping application has an option to load layers from any Geoserver layers. However, as far as I know so far mapfish requires hosts to be written in the config.yaml to be white-listed. This means printing will not work for layers added by a user from another server.

hosts: [*]

Is it possible to have a wildcard to allow all hosts or at least enable to set a variable so that hosts are white-listed dynamically?


texts: Geodjango insert statement: psycopg2.IntegrityError: violates check constraint

texts:
Why is the following INSERT statement throwing this error?

INSERT INTO domes_place (placeid, structidx, oscode) VALUES ('10', '1', 'POINT(51.121282945772 1.3138398501222)');
psycopg2.IntegrityError: new row for relation "domes_place" violates check constraint "enforce_srid_oscode"

This is the table definition in postgres:

           Table "public.domes_place"
   Column   |          Type          | Modifiers
------------+------------------------+-----------
 placeid    | integer                | not null
 structidx  | integer                |
 oscode     | geometry               |
Check constraints:
    "enforce_dims_oscode" CHECK (st_ndims(oscode) = 2)
    "enforce_geotype_oscode" CHECK (geometrytype(oscode) = 'POINT'::text OR oscode IS NULL)
    "enforce_srid_oscode" CHECK (st_srid(oscode) = 4326)

Isn't this POINT in WGS84 - i.e. in SRID 4326?