texts:
Is there a way to clip a linestring based on distance?
For example, linestring from a to b is 150 meters. I would like to get the linestring from a to a point that is is 75 meters away from a following the a-b linestring.

What I have thought until now is to get the line NpointCount and start adding together small line segments until I get close to 75 meters, but still, how am I supposed to add the final line segment?

I hope I am making sense here..

thanks in advance!


texts:
My goal: create viewsheds for observers which also incorporate the heights of buildings which could block observer view. This is using archaeological data, and is representing time change -- so I need to make several versions of the DEM, representing different periods of time when particular buildings would be erect or torn down.

I made raster files representing the extent and height of each building. They're really simple -- just rectangles with a single height. I used the Mosaic to New Raster tool to add three building rasters to a single raster file. Then, I used Mosaic to New Raster again, and added the collective building raster to the DEM (7806x9732) -- when I initially tested this using three buildings and the DEM, it worked fine.

Now, I am running into issues of the tool crashing with non-tested material. I've tried mosaicing each individual building (e.g. DEM + Building 1, DEM and Building 1 + Building 2, etc.) and it crashes. I've remade some building files to check and make sure they aren't corrupt.

I've tried a Con is Null statement as well, and it's crashing:

Con(IsNull("BUILDING"), "DEM", "DEM"+"BUILDING")

I'm getting 999999 errors. Or ArcMap just quits.

Any ideas?

If it's a memory issue, I'm willing to use the polygon version of the extents as well -- all I really need to do is change the values of the DEM within small rectangular areas - adding the height of the building to the DEM.


texts:
I try to download a web map from arcgis online into arcgis desktop. From the site I get a file named item.pkginfo and arcgis desktop opens it. Unfortunately it only gets me the basemap image and not the features.

According to the knowledgebase article 3824 from Esri, one can reach the data via the url after the  tag but what i get is something looking like a long string of json-like data. (The esri-page doesnt say anything about the format one get.)

Haw anyone tried to extract data from arcgis online the same way and can tell the data-format of that string?

The datastring is this: http://www.arcgisonline.com/sharing/content/items/4dde65a4ffd846ff862ed6c508853e37/data?f=pjson for reference.

I admit its a rather vague question but I can't really formulate it better.


texts: Finding the overlap between two rasters based on pixel values

texts:
I'm trying to find the intersection of two rasters.  Raster A has the probability of a fish existing in each pixel and Raster B has the probability of a different fish existing in the same space.  I need to find out where both rasters are showing a probability of 50% or greater, 60% or greater, etc in the same area.  I can see that the probability of 50% for each raster overlaps, but I need to quantify the area.

I don't work with raster often.  I tried converting to vector (polygons) and using the intersect tool, but the resulting polygons were very funky looking and did not represent the total area of overlap.

I'm using ArcGIS 10.6.


texts: QGIS: Buffer around external border of a country based on federal state borders

texts:
based on the question "How to add state polygon information to a city point?" I have an additional question:

With "join attributes by location" I got a layer that constist of the city names with the corresponding ferdale state names. There is just one problem: Lots of the cities are located on the coast and maybe they are not accurate located.

My idea is to create a buffer-zone around the external border.
Is there a possibility in QGIS?


texts: How to use a custom function inside GEE ee.Reducer?

texts:
I have generated a segmentation using SNIC on GEE and calculated some statistics/features for each segment using ee.Reducer, like mean, std, median, etc. I want to calculate more features based on my own equations. However, I didn't figure out how I can use a custom function inside the ee.Reducer. If there is no functionality for doing so, what optimal approach do you recommend me to calculate my custom statistics for each segment?

The script I'm using is something like the following:

var imageCollection = ee.ImageCollection('USDA/NAIP/DOQQ');
var geometry = /* color: #0b4a8b */ee.Geometry.Polygon(
        [[[-121.89511299133301, 38.98496606984683],
          [-121.89511299133301, 38.909335196675435],
          [-121.69358253479004, 38.909335196675435],
          [-121.69358253479004, 38.98496606984683]]], null, false);

var cdl2016 = ee.Image('USDA/NASS/CDL/2016');

var bands = ['R', 'G', 'B', 'N']
var img = imageCollection
    .filterDate('2015-01-01', '2017-01-01')
    .filterBounds(geometry)
    .mosaic()
img = ee.Image(img).clip(geometry).divide(255).select(bands)
Map.centerObject(geometry, 13)
Map.addLayer(img, {gamma: 0.8}, 'RGBN', false)

var seeds = ee.Algorithms.Image.Segmentation.seedGrid(36);

// Run SNIC on the regular square grid.
var snic = ee.Algorithms.Image.Segmentation.SNIC({
  image: img,
  size: 32,
  compactness: 5,
  connectivity: 8,
  neighborhoodSize:256,
  seeds: seeds
}).select(['R_mean', 'G_mean', 'B_mean', 'N_mean', 'clusters'], ['R', 'G', 'B', 'N', 'clusters'])

var clusters = snic.select('clusters')
Map.addLayer(clusters.randomVisualizer(), {}, 'clusters')
Map.addLayer(snic, {bands: ['R', 'G', 'B'], min:0, max:1, gamma: 0.8}, 'means', false)

// Compute per-cluster stdDev.
var stdDev = img.addBands(clusters).reduceConnectedComponents(ee.Reducer.stdDev(), 'clusters', 256)

// Area, Perimeter, Width and Height
var area = ee.Image.pixelArea().addBands(clusters).reduceConnectedComponents(ee.Reducer.sum(), 'clusters', 256)

var minMax = clusters.reduceNeighborhood(ee.Reducer.minMax(), ee.Kernel.square(1));
var perimeterPixels = minMax.select(0).neq(minMax.select(1)).rename('perimeter');
Map.addLayer(perimeterPixels, {min: 0, max: 1}, 'perimeterPixels');

var perimeter = perimeterPixels.addBands(clusters)
    .reduceConnectedComponents(ee.Reducer.sum(), 'clusters', 256);

var sizes = ee.Image.pixelLonLat().addBands(clusters).reduceConnectedComponents(ee.Reducer.minMax(), 'clusters', 256)
var width = sizes.select('longitude_max').subtract(sizes.select('longitude_min')).rename('width')
var height = sizes.select('latitude_max').subtract(sizes.select('latitude_min')).rename('height')


texts: Using Azure Blob storage to save raster and vector data and use it with GeoServer

texts: My client has large datasets store in Azure Blob storage. Is there any way to use them directly via Azure Blob storage? I had a look at Azure GWC plugin but it is just allowing me to cache the existing raster image in blob storage.


texts: Export Google Earth Placemarks to shapefile

texts:
We were provided with a few geotagged photos, from those photos we were expected to mark more classification points around the region for the classes like water, built-up and barren land.

I was able to import those images into google earth using a KML file format using this link. I further went on to mark and label required points around the region. This was done using google earth placemarks. I exported the KML file which has this format.

<Placemark>
    <name>Barren land</name>
    <LookAt>
        <longitude>78.01311981239444</longitude>
        <latitude>13.13658990287736</latitude>
        <altitude>0</altitude>
        <heading>0.0456062423736624</heading>
        <tilt>30.00023885976806</tilt>
        <range>264.4081533793532</range>
        <gx:altitudeMode>relativeToSeaFloor</gx:altitudeMode>
    </LookAt>
    <styleUrl>#msn_track</styleUrl>
    <Point>
        <gx:drawOrder>1</gx:drawOrder>
        <coordinates>78.01311981239444,13.13658990287736,0</coordinates>
    </Point>
</Placemark>

Now I want to create a classifier with these labels.

My query is how to convert Google Earth Placemark KML to shapefile with proper labels for classification using a code in python and not a tool?


texts: Adding polygons with excessively complex coordinates to Google Maps

texts:
I have a big JSON file, from a converted shp file (through http://mapshaper.org/) of region coordinates for Australian Local Government areas.

The problem is that the region areas are super complex and result in lines of coordinates over 200,000 characters long. The JSON file for one state (NSW) is 55mb in size.

These are obviously all too big to be used in any web application.

As I want to draw the polygons on Google Maps, I need to come up with a better way to represent these coordinates. Maybe even a way to simplify the shapes, down to 10% of their points or something?


texts:
I'm trying to plus two raster images through the python interface of QGIS, i want to use the saga raster calculator, but i got one error but i don't know why:
i'm doing the following calculation (the extent of the input date is the same of the raster created by grass in this script

##[test]=group
##test=name
##dem=raster
##output_rast=output raster

GAUSS_1=processing.runalg('grass7:r.surf.gauss', '0.0','15','219818.5071000000461936,245558.5071000000461936,6140993.5039999997243285 ,6177833.5039999997243285',30.0,None)
DEM_Gauss=processing.runalg('saga:rastercalculator', dem,[GAUSS_1],'a+b',False,7,output_rast)

but the code fails, anyone can see the error from this code?


texts: PostgreSQL 12 crash trying to add PostGIS

texts:
When I try to execute CREATE EXTENSION postgis;

I get this error :

geoall=#CREATE EXTENSION postgis;

server closed the connection unexpectedly
        This probably means the server terminated abnormally
        before or while processing the request.
The connection to the server was lost. Attempting reset: Failed.

these are the logs for PostgreSQL

2020-04-15 12:53:23.732 EDT [177498] DETAIL:  Failed process was running: CREATE EXTENSION postgis;
2020-04-15 12:53:23.733 EDT [177498] LOG:  terminating any other active server processes
2020-04-15 12:53:23.734 EDT [177787] db@geoall WARNING:  terminating connection because of crash of another server process
2020-04-15 12:53:23.734 EDT [177787] db@geoall DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, becauther server process exited abnormally and possibly corrupted shared memory.
2020-04-15 12:53:23.734 EDT [177787] db@geoall HINT:  In a moment you should be able to reconnect to the database and repeat your command.
2020-04-15 12:53:23.734 EDT [177786] db@geoall WARNING:  terminating connection because of crash of another server process
2020-04-15 12:53:23.734 EDT [177786] db@geoall DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, becauther server process exited abnormally and possibly corrupted shared memory.
2020-04-15 12:53:23.734 EDT [177786] db@geoall HINT:  In a moment you should be able to reconnect to the database and repeat your command.
2020-04-15 12:53:23.737 EDT [177783] WARNING:  terminating connection because of crash of another server process
2020-04-15 12:53:23.737 EDT [177783] DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, because anotherr process exited abnormally and possibly corrupted shared memory.
2020-04-15 12:53:23.737 EDT [177783] HINT:  In a moment you should be able to reconnect to the database and repeat your command.
2020-04-15 12:53:23.758 EDT [177498] LOG:  all server processes terminated; reinitializing
2020-04-15 12:53:23.777 EDT [179200] LOG:  database system was interrupted; last known up at 2020-04-15 12:52:56 EDT
2020-04-15 12:53:23.777 EDT [179201] postgres@geoall FATAL:  the database system is in recovery mode
2020-04-15 12:53:23.823 EDT [179200] LOG:  database system was not properly shut down; automatic recovery in progress
2020-04-15 12:53:23.828 EDT [179200] LOG:  redo starts at 2/996F5960
2020-04-15 12:53:23.885 EDT [179200] LOG:  redo done at 2/99945FC0
2020-04-15 12:53:23.914 EDT [177498] LOG:  database system is ready to accept connections

and I can't install PostGIS.

Version of Linux:

Linux  5.4.0-4-amd64 #1 SMP Debian 5.4.19-1 (2020-02-13) x86_64 GNU/Linux


texts:
Are there any online statistics for the popularity/usage of web mapping APIs? A diagram would be sufficient, if it can be published in a paper. I have only found one source: https://trends.builtwith.com/mapping.
The problem with this source is, the license states that it cannot be published.

You may view, download for caching purposes only, and print pages from the website, provided that:
(a) you must not republish material from this website (including republication on another website), or reproduce or store material from this website in any public or private electronic retrieval system;
(b) you must not reproduce, duplicate, copy, sell, resell, visit, or otherwise exploit our website or material on our website for a commercial purpose, without our express written consent;

Are there any citable sources out there for web mapping trends? Or did I misunderstand the conditions, and it can be published with proper citing?


texts: Extracting By Mask list of raster files using ArcPy?

texts:
I'm trying to mask a group of raster files and save the new rasters in a different folder using the same name of the original files.

It doesn't sound as a difficult task, but I'm having trouble saving-naming each masked raster inside the "for" loop. Here is what I have so far:

import arcpy
from arcpy import env
from arcpy.sa import *

env.workspace = "C:/rasters/threshold"
mask = "C:/GIS/mesoamerica.tif"

rasterlist = arcpy.ListDatasets("*", "Raster")
    for i in rasterlist:
        outExtractByMask = ExtractByMask(i, mask)
        outExtractByMask.save("C:/SIG/MelasCA_30runs_avg/threshold/mesoamerica/" i)

I need to find a way to use the "i" variable in the outExtractByMask.save. Any idea?


texts: Georeferencer GDAL Plugin Not Present in Manage/Install Plugins

texts: I have a student who has installed QGIS 3.x on her Mac OS and cannot locate the Georeferencer GDAL plugin. It does not appear in "installed" (as it should) nor via search in the plugins manager. An uninstall/reinstall of QGIS has not made it appear. Any tips?


texts:
I have collected the signal strengths for the mobile operator at some places in my college campus. Now I want to spatially interpolate this data for the entire campus. I am using PostGIS.

I am planning to use the Inverse Distance Weighting (IDW) method for the interpolation. Can someone guide me how can I employ this by using PostGIS. Are there any inbuilt functions available or will I have to write the separate algorithm for this.


texts: Creating plots for each hour of each day with Python

texts:
I am trying to create an animated map (by generating multiple plots) of road traffic throughout a week, where the thickness of roads is represented by the volume of traffic at a specific time of day.
This is sort of what I'm looking for (but for each hour of each day):

The data has a structure that looks like this:
HMGNS_LNK_ID                 geometry    DOW  Hour   Normalised Value
         2     MULTILINESTRING ((251...    1     0   0.233623
         2     MULTILINESTRING ((251...    1     1   0.136391
         2     MULTILINESTRING ((251...    1     2   0.108916

DOW stands for 'day of the week' (1 = Monday) and so for every Hour of each of the 7 days I want to plot the map with roads' thickness by the value Normalised Value.
I encounter a problem that when trying to loop with this code:
for dow in df['DOW']:
    fig, ax = plt.subplots(1)
    day_df = df[df['DOW']==dow]

    for hour in day_df['Hour']:
        day_hour_df = day_df[day_df['Hour']==hour]
        day_hour_df.plot(ax=ax, linewidth=day_hour_df['Normalised Value'])
        plt.savefig("day{}_hour{}.png".format(dow, hour), dpi = 200, facecolor='#333333')

The problem is that the figures are saved only for day 1, so until day1_hour_23 and after that, it comes back to day1_hour0 and overwrites the plot with something new. I can't figure out why it stops at DOW 2.
I'm not even sure if the data structure is correct.
The full code is in my GitHub repository.


texts: Point locator on polygon for onclick of elevation profile graph

texts:
I am trying to display a point locator on polygon path, when onclick of elevation profile graph.I tried the following code with google api in openlayers library, the graph is showing but pointer is not displaying. Please suggest me how to add a pointer on polygon path.

var elevator = new google.maps.ElevationService;
            displayPathElevation(path, elevator, map);

function displayPathElevation(path, elevator, map) {

    // Display a polyline of the elevation path.
   /* new google.maps.Polyline({
      path: path,
      strokeColor: '#0000CC',
      strokeOpacity: 0.4,
      map: map
    });*/

    // Create a PathElevationRequest object using this array.
    // Ask for 256 samples along that path.
    // Initiate the path request.
    elevator.getElevationAlongPath({
      'path': path,
      'samples': 150
    }, plotElevation);
  }

  function plotElevation(elevations, status) {
    var chartDiv = document.getElementById('elevation_chart');
    if (status !== 'OK') {
      // Show the error code inside the chartDiv.
      chartDiv.innerHTML = 'Cannot show elevation: request failed because ' +
          status;
      return;
    }
    // Create a new chart in the elevation_chart DIV.
    var chart = new google.visualization.ColumnChart(chartDiv);

    // Extract the data from which to populate the chart.
    // Because the samples are equidistant, the 'Sample'
    // column here does double duty as distance along the
    // X axis.
    var data = new google.visualization.DataTable();
    data.addColumn('string', 'Sample');
    data.addColumn('number', 'Elevation');
    for (var i = 0; i < elevations.length; i++) {
        displayElevation.push(elevations[i].elevation)
      data.addRow(['', elevations[i].elevation]);
    }

    // Draw the chart using the data within its DIV.
    chart.draw(data, {
      height: 150,
      legend: 'none',
      titleY: 'Elevation (m)'
    });

}


texts: How do I download all features with a particular tag?

texts:
I'm trying to get all the different features(point, line & polygon) which are tagged with locality, from OSM data from a metropolitan sized area

I'm presently working with minimal computer resources so I can't upload all osm data from this region in a PostGIS database.

Is there any online website or service which allows the download of OSM data after filtering it by tags?


texts: Exporting several maps with excel sheet using Arcpy

texts:
I try to export 20 maps with this python code in ArcGIS 10.3.1 for Desktop. Each map has excel sheet (not attribute table) in the layout view. The excel sheet file located in G:\Users\test\Downloads and imported from there:

import arcpy,os
from arcpy import env

env.workspace = r"G:\desktop\Project"
for num, mxd in enumerate((arcpy.ListFiles("*.mxd")),start = 1):
    print num,mxd
    mapdoc = arcpy.mapping.MapDocument(r"G:\desktop\Project\\" + mxd)
    arcpy.mapping.ExportToJPEG(mapdoc,
    r"G:\desktop\Project\ou\\" + os.path.basename(mapdoc.filePath).replace('.mxd','') +".jpg",resolution=196)
mapdoc.save()
del mxd

This code works fine, but when i insert object- excel sheet with those two options:

when i use the python code i get only white rectangle with no data table:

i red this , but didn't found any solution:

Attribute table added to layout is missing when mxd is exported to jpeg using ArcPy?
Arcgis table and Excel sheet dynamic linking


texts: Am trying to edit a line shapefile using a spline tool in QGIS 2.4.0 .Its drawing but the ok button is inactive. what could be the problem? Generally am trying to create a curved line.


texts: Is there an optimal algorithm to find the longest shortest path in a network?

texts:
I have a large set of linear networks and I would like to find the two ends of each network that are the most distant from each other along the network (on the image below, it would be D to K). The brute force solution to this problem is to compute the shortest path along the network for each pair of origin, but I have hundreds of networks with thousands of ends, so computing each possible path is quite heavy.

Is there an optimal way to compute this without using the brute force ? Can I exclude some points based on some clever rules ?

EDIT: I've added an illustration of the longest path mentioned by @Alex Tereshenkov in order to clarify my question.  The black path is the result of the longest path algorithm (longest path without repeating any vertices). In my case, imagine that you enter the network from any of the letters and you need to drive to another letter as fast as you can. Which two letters are the most difficult to join ?


texts: Mapping carbon sequestration potential as green infrastructure function?

texts:
I have a small project where I have to do mapping of potential areas for carbon sequestration in the Adelaide Metropolitan area in South Australia. I have no idea what data do I need.

So far looking at literature I understand I can do it through NDVI but have no further idea.

My project is simple, just showing the potential carbon storage area only, not to calculate C02.

Can anyone help me in finding methodology or relevant literature.


texts: Directions based on polygon centroid and average/median point of adjacent polygon's touching border

texts:
I have a layer containing multiple polygons, and I am trying to find a way to identify the directions of the adjacent polygons of a certain polygon. The directions would be based on the main polygon centroid (red dot) and the average/center point (blue dot) of the touching polygon's line:

Is there any script, plugin, or algorithm that accomplishes this, preferably one that can do it in batch?


texts: Using twitter coordinates in OSGeo library

texts:
I have twitter data with coordinates such as: [-7.2377761, 112.73091137]

I'm working with a python geo library called osgeo. The sample code indicates creating a point like this.

point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
point.SetPoint(0, 474595, 4429281)

The values above do not look like lat/lon. How do I perform the conversion?


texts: How to find and quantify duplicate vertices on a polygon/polyline in ArcGIS 10.0

texts:
I need to clean up topological errors in various shapefiles and databases. I did read that the Repair Geometry tool could clean duplicate vertices but how can I know the number of errors before and after using this tool, to ensure these errors were repaired?

I'm using ArcGIS 10.0. I have no programming skils.


texts: arcPy findAndReplaceWorkspacePath Won't Update Rasters

texts:
I'm using ArcGIS 10.3. I have an MXD that contains symbolized global raster images. I have a GP tool that masks these rasters to a country or state and exports the images. I'm trying to update the source background layers and rasters each time a new country is chosen. The layer names are all generic, but the source files all have the country or state name tagged on.  findAndReplaceWorkspacePath works just fine with the background layers (shapefiles), but it won't update the rasters. I've confirmed the old and new sources with print statements, it's just not changing the source. The only time I get images is if I use the default state (Texas) set in the MXD. If I choose any other location, I just get a blank country or state outline image.

Here is my code snippet:

 layer_list = ["Precip","High_Temp","Temp_Departure"]
 for filename in layer_list:
   layer_name = filename
     for lyr in arcpy.mapping.ListLayers(mxd):
       if lyr.supports("DATASOURCE"):

          if layer_name == str(lyr.name):
            old_source = os.path.dirname(lyr.dataSource)

            new_source = country_temp_path

   lyr.findAndReplaceWorkspacePath(old_source, new_source, False)

return()


texts: Updating Feature Service that ArcGIS Online Feature Layer points at

texts:
I have a Feature Service hosted on my ArcGIS Server and exposed to the internet (something like this Layer: Rainfall(NOAA) (ID: 0)). I added this Feature Service to my AGOL org by going to Content >  Create > Feature Layer > From URL and entering the REST URL of my Feature Service (I ended up with something like this: https://geodecisions.maps.arcgis.com/home/item.html?id=caa1b120e1ac4fa497866e0ba4cdcd9e). Now I would like to change the Feature Service to a different service on my ArcGIS Server while keeping the same Feature Layer item id. Thie would allow me not to have to update the Feature Layer URL in a few applications that I have that use this Feature Layer.

Is this possible?

I have looked around in the Feature Layer settings and tried the ArcGIS Online Assistant but haven't found anything that works.


texts: CRS options don't appear on Macintosh QGIS version 3.10

texts:
I am teaching a basic GIS class right now, having my students use QGIS on their own computers.
Three of my students, running QGIS 3.10 (on Macs with OS 10.15), are having a problem that there are no options for selecting a CRS either for a project or layer.
When they click on project properties (or layer properties or layer export window), there are simply no loaded CRS options in the window.


texts: Reconciling vertex modification conflicts in Geodatabase Versioning?

texts:
Is it possible to edit the same feature in ArcSDE in two different versions

... and reconcile/post version A's edits to the parent,

... and then reconcile/post version B's edits of same feature,

... and then see which vertex modification conflicts arise,

... and keep all vertex modifications that do not conflict?


texts: Why does gdal2tiles sometimes produce a tileset where the image does not fill the tiles?

texts:
I'm using gdal2tiles on a vrt that I've warped from Mercator to QSC. I'm using the 'raster' profile. The resulting tiles look correct but there transparent strips on the upper and right edges of the right-most and top-most tiles:

The above image shows tiles 1,1,0 and 1,1,1

My call to gdal2tiles looks like this:

gdal2tiles.py --profile=raster -z 0-4 face.tif outFilePath


texts: Allow to insert only one unsaved marker on Google Map

texts:
I use a google Map to collect point data, that is saved to a Postgis database.
With dblclick the user adds a marker, that is draggable and has a listener that opens (onclick) a modal dialog window to enter additional data. When the data is saved to the database (via php-script), the standard marker is replaced by a marker that reflects the type of point acquired. The user can than add the next  point.

All this works fine, I only have the problem, that the dblclick-event can be fired a few times, so that the user gets confused with several new markers.

The addListenerOnce doesn't help me since the user usually needs to add a bunch of new points. Also this old post didn't help googlemaps-v3-api-create-only-1-marker-on-click since it is intended to only use ONE point.

How can I prevent the dblclick-event to add new markers, as long, as the first one isn't saved?

   var marker;
   var modalOnClick;

   //Add listener to map, to insert new Marker with modal
   google.maps.event.addListener(map, "dblclick", function(event){
     marker = new MarkerWithLabel ({
      position: event.latLng,
      map: map,
      icon:mmcross,
      draggable: true,
        labelContent: '',
        labelAnchor: new google.maps.Point(15,33),
        labelClass: "labels", // the CSS class for the label
        labelInBackground: false
     });

    map.setCenter(event.latLng);
      //Open Modal Form to insert data
    modalOnClick = google.maps.event.addListener(marker, 'click', function() {
      $('#myModal').modal('show');
    });
   });

// Function to save data
function saveData(){
  ...

    $('#myModal').modal('hide');
    google.maps.event.removeListener(modalOnClick);
}; // END saveData()


texts: Error in sde.next_rowid values through PostgreSQL Geodatabase

texts:
I have a spatial data model implemented in PostgreSQL / PostGIS with several alphanumeric and spatial tables. The spatial data is edit mainly through QGIS.

As I want to be able to edit the spatial data in ArcGIS environment and publish a Feature Service in ArcGIS Server, I am trying to implement a solution where I have a column as primary key (PK) managed by PostgreSQL and an OBJECTID column to function as PK in the ArcGIS environment managed by SDE (as solution proposed in: Feature Service from plain PostgreSQL/PostGIS in ArcGIS Server 10.2)

To be able to still edit the spatial data through QGIS I had to implement a trigger function that will fill the OBJECTID (SDE PK) field with the appropriate value. In order to do this, I am using the sde.next_rowid function, but the values returned by this function are different (much higher) than the ones attributed to new features created in ArcGIS environment, as if they were using different sequence numbers.

I know that sde.next_rowid function is using the ids defined in the [table_owner].i[sde_id_from_table_registry], but does anyone have idea where SDE is getting the ids for new features created in ArcGIS environment?

Example:

When I run the this SQL statement in PostgreSQL:

SELECT sde.next_rowid('bdsiguser', 'pg_npolicia_pt');

I get a value of 144150.

But, when I edit the same feature class (bdsiguser.pg_npolicia_pt) in ArcGIS and create a new feature, it gets the objectid 143754!

And the next oid values, through each of the methods would be sequencial, i.e., 144151 and 143755 accordingly.

Note: I am using ArcGIS 10.2.2 and PostgreSQL 9.2.2 with PostGIS 2.0.6 (ESRI's recommended PostgreSQL/PostGIS version for this ArcGIS version)


texts: Drawing a map using time as metric instead of distance in QGIS

texts:
I am used to drawing maps in QGIS and overlaying on top them a road network map showing the distance between places, its a useful way of showing information about roads.

But I am aware that its not showing everything,  some classes of road  allow much faster travel e.g. 10 minutes driving on a single track  country road and I might cover two miles, were as on a motorway I might cover 10 miles.

Is there anyway I can use this information do distort either a shapefile or raster image of a generic to map to show the faster areas compared to the slower areas?  I have a set of X/Y coordinates and the time to travel between them.

I am using QGIS to display my data.  I normally keep my data in a PostGIS database.


texts: Looking for documentation on GRASS vectors/topology

texts:
I am looking for documentation (including pictures) of the rationale behind GRASS topology that answers questions like: "What is a topology according to GRASS? How are topologies inferred from geometries, how are nodes, segments and vertices defined?"

So far the internet has let me down in the search, but maybe I missed something.
The best reference I have found is: https://grass.osgeo.org/grass72/manuals/vectorintro.html
but it doesn't really dive into the working of the topology. Ideally there should be something like: https://bost.ocks.org/mike/topology/ (which does not work the same as grass topology).


texts: Is there a way to convert 2D features to "3D" features in QGIS?  I know there is a way to do this in ArcMap by using an attribute value in the attribute table, but I am looking to do this with shapefiles in QGIS.


texts:
I'm working with an area spanning two UTM zones. The magazine I will submit my work to needs the projection in UTM.

Is it possible to combine two neighbouring UTM zones to present a map/data that spans them?


texts: Creating multiple shapefiles from one date within same for loop using ArcPy?

texts:
I have imported North American wildfire incident data a comma delimited text file containing one line for each fire incident. Each fire incident has a latitude, longitude coordinate pair separated by commas along with a confidence value.

My first task was to create a shapefile, and make the data usable to import into the attribute table of the shapefile I created called Incidents.

My new task is to repeat this process in a new shapefile within my MDB. The shapefile is called NewFires by user input. My issue is that when I run the code a second time, I get the error:

  Traceback (most recent call last):   File "C:\Program Files (x86)\Wing
  IDE 101 4.1\src\debug\tserver_sandbox.py", line 29, in
  File "C:\Program Files
  (x86)\ArcGIS\Desktop10.3\arcpy\arcpy\management.py", line 1807, in
  CreateFeatureclass
      raise e arcgisscripting.ExecuteError: ERROR 000354: The name contains invalid characters Failed to execute (CreateFeatureclass).

My second shapefile starts with 'if 'Confidence' in fire:'

import arcpy

work = raw_input("Enter the full path of WildlandFires.mdb: ") # C:\Scripts\Lab 7 Data\WildlandFires.mdb
arcpy.env.workspace = work # Set the workspace to the geodatabase
arcpy.env.overwriteOutput = True

iFile = raw_input("Enter the full path of wildfire text file: ")# input text file C:\Scripts\Lab 7 Data\NorthAmericaWildfires_2007275.txt
output = raw_input("Enter the name of the output feature class: ")

f = open(iFile, 'r')# input text file C:\Scripts\Lab 7 Data\NorthAmericaWildfires_2007275.txt in "read" mode

lstFires = f.readlines() # Read the lines of fire data from the input file
f.close() # close the file

fields = ["SHAPE@", "CONFIDENCE"] # represents the field names for each row in the feature class

# Process: Create Feature Class
arcpy.CreateFeatureclass_management(work, "Incidents.shp", "POINT", "C:\\Scripts\\Lab 7 Data\\WildlandFires.mdb\\FireIncidents")

cursor = arcpy.da.InsertCursor("Incidents.shp", fields) # add the new points that you will create to the "FireIncidents" feature class

cntr = 0 # counter variable

for fire in lstFires:
        if 'Latitude' in fire: # Skip the header
            continue
        pnt = arcpy.Point() # Create a new Point object
        lstValues = fire.split(',') # create a list of values "lstValues" based on the type of delimiter
        str_int = lstValues [-3:]
        str_int = list(map(float, str_int)) #turns string to integers
        latitude = str_int[0] # Latitude
        longitude = str_int[1] # Longitude
        confid = str_int[2] # Confidence Value

        pnt.X = longitude # Assign the X and Y properties of the Point object
        pnt.Y = latitude

       if 'Confidence' in fire:  # Skip the header
            continue
        output = arcpy.Point() # Create a new Point Object
        arcpy.CreateFeatureclass_management(work, "NewFires", "POINT", "'C:\\Scripts\\Lab 7 Data\\Incidents.shp'")
        arcpy.AddField_management(NewFires_shp, "CONFIVALUE", "FLOAT")
        field = ["SHAPE@", "CONFIDENCE"] # represents the field names for each row in the feature class
        row = [pnt, confid] # Create a new row for the feature class
        cursor.insertRow(row) # Insert the new row to the feature class

        cntr = cntr + 1 # update the counter
        print "Record # " + str(cntr) + " written to feature class"

del cursor # release the cursor lock on the feature class


texts: What gestures make sense in a GIS application?

texts: Have been testing Microsoft Beta APIs for kinect with the kinect Toolkit. To create an appliation that catch a users gesture is very easy. As a GIS developer it feels like it would integrate perfectly with ESRI ArcGIS for WPF, to build for example kiosk applications where a user can interact with a map. What gestures would make sense for a user to navigate a map?


texts:
I try to merge two DEM data sets with python, using gdal_merge, but it doesn't work. I install like this:

from osgeo import gdal

But this doesn't work:

import gdal_merge as gm

I get the error message:

  No module named 'gdal_merge'

I am entirely new to Python


texts: Identifying the closest polygon to a point dataset

texts:
I have two large data sets one point and one polygon.

What I was looking to do was identify which polygon is closest to each point.

Hopefully I was looking for a automated process as the point data set is very large. Ideally I would like a column in the attribute table of the point data set identifying the name/reference of the nearest polygon.


texts: Replacing shapefile with feature class in geodatabase using ArcPy?

texts:
I have a shapefile which I turned to feature class to perform a rename on some fields.

arcpy.FeatureClassToGeodatabase_conversion(r'C:\Users\user\Desktop\Dty\Dty.shp', r'C:\Users\user\Desktop\gw\gb.gdb' )

#SET ENVIRONMENT TO THE GDB
arcpy.env.workspace = r'C:\Users\user\Desktop\gw\gb.gdb'
#rename
arcpy.AlterField_management('Dty','KYR_DATE1','KYR_DATE')

#Now turn this back and replace the initial shapefile
arcpy.env.overwriteOutput = True
arcpy.FeatureClassToShapefile_conversion('Dty',r'C:\Users\user\Desktop\Dty')

last line doesn't replace the original shapefile, it creates Dty_1.shp
 What is the correct way?


texts: Raster displays on different locations on windows and ubuntu

texts:
I have worked on some continental scale raster and vector datasets and all rasters are saved as Netcdf´s.

This was done by converting rasters with different timestamps and variables to a Netcdf file.

When working on Ubuntu (QGIS 2.14) all vectors and rasters are in the right place but when loading the same data in windows (QGIS 2.16) the rasters get flipped and move out of place.

I have been using EPSG:4326.
Below follows and example of the metadata for a file with one raster layer (Which also displays wrongly).

Driver
GDAL provider
HDF5Image
HDF5 Dataset
Dataset Description
HDF5:"G:/Aquaspace/Final/Raster/BATH/BATH3.nc"://BATH
Conventions=CF-1.5
GDAL=GDAL 1.11.3, released 2015/09/16
GDAL_AREA_OR_POINT=Area
history=Fri Aug 26 17:48:12 2016: ncatted -O -a long_name,BATH,o,c,bathymetry BATH.nc Fri Aug 26 17:47:44 2016: ncrename -v Band1,BATH BATH.nc Fri Aug 26 17:46:35 2016: GDAL CreateCopy( BATH.nc, ... )
Band 1
BATH_grid_mapping=crs
BATH_long_name=bathymetry
BATH__FillValue=-32767
STATISTICS_MAXIMUM=3800
STATISTICS_MEAN=-28506.247088649
STATISTICS_MINIMUM=-32767
STATISTICS_STDDEV=10787.572663429
Dimensions
X: 21601 Y: 10801 Bands: 1
No Data Value
*NoDataValue not set*
Data Type
Float32 - Thirty two bit floating point
Pyramid overviews
Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
0.0000000000000000,-10801.0000000000000000 : 21601.0000000000000000,0.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet

Does anyone know what could cause this problem?


texts: Most efficient way to select a given layer programmatically in QGIS

texts:
this a question about code efficiency, not asking how to loop through layers. I'm getting the 'this question appears subjective...' warning but there must be a tried and tested method by someone more adept at PyQGIS programming than me! I don't intend to start a discussion, just hopefully someone more experienced can objectively say "do it like this instead" or "that's the best/only way".

I'm only just starting to use QGIS and PyQGIS to develop...stuff. Most of the operations I will be doing, and indeed all of them up until this point, will be using some layer that's in the map registry. So every function I have so far starts like this:

def some_function(layer, parameter, another_parameter):
    for lyr in iface.mapCanvas().layers():
        if lyr.name() == layer:
        # and then whatever it is I need to do or get from this layer

Is this the most efficient way of going about this, or is there another way of proceeding without looping through all layers and then having the if statement to ID the layer. Clearly the more layers you have open, the more computationally expensive this would be, perhaps there is a way of selecting the layer by its name directly?


texts: Drag points on lines with snapping in QGIS 2.18

texts:
I have a point layer and a line layer coming from PostGIS. Points represent sensors (two for each side of the road), lines are GPS tracks. I need to move the sensor locations towards the lines. As I cannot take the closest track in each case, I decided to do it manually because I only have around 100 sensors.

That's easy to handle in QGIS, right? But to my surprise, I'm struggeling to activate snapping when moving points. I've enabled snapping for the corresponding layers and I can see that it would work when creating new features but not for moving them. According to this source, I should use the node tool. While it's possible to snap two points with each others they still not land on the track. Are there any other options I missed?


texts: QGIS 2.8.1 (Windows 64bit). Setting snapping to the "decoration" grid?

texts:
I was wondering if anyone knew of a way to snap digitizing to the grid you can create in View>Decorations>Grid?

I can't see anything in the manual or the menus so I am thinking there probably isn't but it would be really helpful for me if you could so I thought I would ask. I'm not great with python but I'm willing to give it a go if that's what it takes.


texts: Using Arcpy to clip Features in Workspace

texts:
I'm trying to write a program that clips all the features in a workspace.  I've read through several similar threads on this website, but still am having trouble understanding what is wrong with my program.

# This program clips all roads, state roads, railroads, steams, rivers, and
# land use using the Meridian Township boundary.  Each clipped dataset is
# renamed to the base name plus "_Meridian.shp"

import arcpy

arcpy.env.workspace = "F:/data"

clipFeature = 'Meridian.shp'

appendstring = '_Meridian'

featureClassList = arcpy.ListFeatureClasses()
# this gets the base name variable
clipFeatureName = arcpy.Describe(clipFeature)
clipName = str(clipFeatureName.baseName)

try:
    for featureClass in featureClassList:
        featureClassDescription = arcpy.Describe(featureClass)
        featureClassList = arcpy.ListFeatureClasses()
        #this gets the base name variable
        nameString = str(featureClassDescription.baseName)
        arcpy.Clip_analysis(featureClassList, clipFeature, "F:/output")

        # this statement prevents the clip feature from clipping itself.
        if nameString == clipName:
            pass
        else:
            namething = nameString + appendstring
            print namething

except:
    print arcpy.GetMessages(2)

Update: My program does not produce an error message, but when I look in my output folder in arccatalog, it is empty.


texts: Resample rasters in for loop using arcpy/python

texts:
I have a folder of annual raster land cover maps that are 56x56m spatial resolution. I would like to resample them to 30x30m using a for loop in that folder. I think I am making a small error but unsure of what it may be. The function works but is only able to resample and save one file then it stops.

import arcpy, os
from arcpy import env
arcpy.env.overwriteOutput = True

env.workspace = r"C:\Users\Documents\rasters"
outFolder = r"C:\Users\Documents\new_rasters"

rasterList = arcpy.ListDatasets("*","Raster")
for inRaster in rasterList:
    outName = os.path.join(outFolder,str(inRaster))
    outResampled = arcpy.Resample_management(inRaster, outName, "30",  "NEAREST")
    outResampled.save(outName)


texts: Error getting text attributes from shapefile with python / pyshp

texts:
I have a shapefile, containing only points in XY-coordinates and a string attribute for each point. I would like to read the coordinates and the attribute for postprocessing in python, but can't get it to work.

I have tried the following.

import shapefile
r = shapefile.Reader("subregion_59117")
print r.fields
shapes = r.shapes()
print len(shapes)
print shapes[0].points[0]

gives

[('DeletionFlag', 'C', 1, 0), ['LON', 'N', 16, 8], ['LAT', 'N', 16, 8],
['SUBREGION', 'C', 25, 47]]

59117

[176.0, 63.0]

So I can read the point coordinates, but

records = r.records()
sr = r.shapeRecords()
print len(records), len(sr)

gives just

0 0

meaning, these lists don't contain any elements. So I don't know how to get my hands on the string attributes of the points?

However, if I open the shapefile in QGIS, and save as geojson (a JSON textfile, thus human readable), it contains the features as a list of elements such as:

{ "type": "Feature", "id": 0, "properties": { "LON": 176.0, "LAT": 63.0, "SUBREGION": "Eastern Europe" }, "geometry": { "type": "Point", "coordinates": [ 176.0, 63.0 ] } }

So apparently my shapefile contains all the information, I just don't know how to read it with pyshp.

(So yes, I could always convert my shapefile to geojson with QGIS, and postprocess that with python, but this feels overtly complicated.)

My shapefile is here.

Edit: My shapefile was generated by Saga GIS, and as explained by Tuure Laurinolli's answer, Saga GIS seems to have a bug so the attributes in my shapefile were marked with an invalid flag, and also pyshp seems to be very picky in this, interpreting everything but the correct flag as "deleted" and refusing to read it.

Converting the broken shapefile to shapefile by e.g.

 ogr2ogr file_new.shp file_broken.shp

fixes the problem.


texts: Different clipping using gdalwarp on Python and GIS

texts:
I'm doing some clipping tests with many polygons on different tiles and I noticed that the gdalwarp command returns two different clippings with Python and QGIS, could you tell me why?

import gdal  # version 3.0.4
import os

polygons = #path shp list
rasterfiles = #path .tif list

i=1
for polygon in polygons:
    subdirshp,shp=os.path.split(str(polygon))
    shp= os.path.splitext(shp)

    for raster_path in rasterfiles:
        subdir,file=os.path.split(str(raster_path))
        file= os.path.splitext(file)
        output_path = outputPath + "\\ID_{}_TILE_{}".format(shp[0],str(i)) + ".tif"

        result = gdal.Warp(output_path, raster_path, cutlineDSName=polygon, cropToCutline=True, dstNodata = 0.0)
        i=i+1
        outRaster = gdal.Warp(srcDSOrSrcDSTab=rasterfiles, destNameOrDestDS=output_path +"\\ID_{}_TILE_{}".format(shp[0],str(i)))

        outRaster= None

Command line on sample QGIS:

gdalwarp -of GTiff -cutline ./input_shp.shp -crop_to_cutline -dstnodata 0.0 ./input_raster.tif ./otput.tif


texts: Is it possible to create an ArcGIS library?

texts: Can I create my own ArcGIS library that contains commonly used features inside my own functions/methods? I would like to make a dll that can be reference in my projects for using common features such as Geo-Coding addresses or creating a FieldsCollection or an ITable. Will this be a problem due to licensing and having to verify it?


texts: I am trying to find a program or application that would allow co-workers to view maps and various layers produced with QGIS but not have the ability to alter or modify the data.  It would be nice to produce something more interactive than just sharing a PDF map, any suggestions of a way to proceed??


texts: Is it possible in QGIS print composer to include a numeric scale in an html box?

texts:
For work I always have to create a table with my name, project name, date etc. into my map layouts. I also need to include the scale like "1:150000" in that table.

So far I am solving that by including an html table with the html box function of the print composer. Is there a way to automatically include the right scale? So far I only managed to automatically include the date by using a little javascript. Maybe something similar is also possible for the scale?


texts: How to develop a Toolbox intended to be distributed in multiple version of ArcGIS?

texts:
I am currently developing a Toolbox for users who will use it in different versions of ArcGIS (10+ versions). Is there a recommended way to proceed in that situation? My main concerns about that are:

Someone had to face that issue before me, how did they overcome the "one version per machine" limitation (ESRI FAQ) ? You can install another version on a virtual machine, but it makes the whole thing more laborious.
How to deal with arcpy version differences in the code? Is there a clean way to handle multiple arcpy versions in the same application? All that I can think of is either split completely the tools for the different versions or design some sort of facade/strategy design pattern.

I'm quite surprised I couldn't find more information on how people deal with that.

For people passing by: the Toolbox must be saved in the oldest version (ArcGIS Help 10.1)


texts:
I have been looking for solutions on the Internet for the last two days without finding answers. Currently, I have a Java web application that displays a map using OpenLayers. The problem is that the map is static, which means that it is entirely build from our images depending on the area we want to display, causing some memory trouble.

We would like to do that dynamically (i.e. send a request to load an image correspondinf to the part of the map we want to display).

I have read many things about Tile Map Service, Web Map Tiling Service, Web Map Server, TileCache, etc... I just don't understand it exactly works. It seems that WMTS seems to be a good choice because of it simplicity and speed.

My questions are, how do I put all this together and what should I use ? What is the architecture I need to set on server side (using OpenLayers on client side if possible) ?

Thank you in advance.


texts: Having trouble selecting polygon from existing layer (using OpenLayers.Control.SelectFeature)

texts:
I am hoping someone can point me in the right direction. I am not a programmer, so that is my first roadblock. ha

The problem:
I have an existing polygon layer and I am trying to select one of the polygons. This has got to be one of the simplest tasks, but I am unable to figure out how to do it. I am trying to follow the methods in the OpenLayers Beginner's Guide, but so far nothing has worked out:

var select_feature_control = new OpenLayers.Control.SelectFeature(
     legal_land_descriptor,
     {
    multiple: false,
    toggle: true,
    multipleKey: 'shiftKey'
     }
);

map.addControl(select_feature_control);
select_feature_control.activate();

Should this not allow me to select a polygon?  I did manage to query the polygons by clicking on them, and now I am just trying to select one of the polygons, Argh!

Any ideas as to what I am doing wrong here?

Thanks!!


texts: Comparison of GDB Table with a database table Comparison

texts:
I have to compare 2 tables one from gdb and another from database. I do not want myself to do a comparison one by one. Do anyone of you have a better idea for doing this?

UPDATE: I have read the documentation. I have tried with the Table Compare. But the output came like this for the below scenario#
Scenario 1) DB table does not have any data.

Has_error, Identifier, Message, Base_value, Test_value, ObjectID

"true", "Table", "Table row counts are different", "1734", "578", -1
"false", "Table", "Tables have same number of fields", " ", " ", -1

Scenario 2) DB table has data 1 time insert.

Has_error, Identifier, Message, Base_value, Test_value, ObjectID

"true", "Table", "Table row counts are different", "578", "2312", -1
"true", "Field", "TestField lengths are different", "8", "36", -1

"false", "Table", "Tables have same number of fields", " ", " ", -1
"false", "Table", "Rows are the same", " ", " ", -1

The field is a date field which I suppose differnt in GDB and Database thats why the field is updated in the file as I think. I may be wrong. But not sure what is the way to use this Table Compare in such cases.


texts: Standard saving .mxd for earlier versions?

texts:
I use two computers, one with ArcGIS 10.2, the other one uses 10.5. I know, i can use "save a copy" to save he project for the older version of ArcGIS, but then, it is just a copy.

Is it possible to use a project in 10.5 and save it for the 10.2 version in a standard way?

Like just clicking on "save" and it always overwrites the 10.2 version although using 10.5?


texts:
I want to do an overlay analysis to find most suitable locations to live in an area.

My criteria are road network, water lines, hospital locations, and schools. My steps are:

Convert all layers to raster using Rasterize
Apply Proximity (raster distance) tool
Apply Ordered Weighted Averaging tool to get the final map.

But it does not display the result due to the error: "The above files could not be opened, which probably indicates that they were not correctly produced by the executed algorithm." I tried this many times but result was same.

I want to know what the wrong with my method or if there is another Vector or Raster based method to perform this type of analysis using QGIS. I am using QGIS 2.2 Valmiera.


texts: Classifying DEM using Height Differences in ArcGIS ModelBuilder?

texts:
I would like to classify DEM into some classes, into a shp file with equal intervals of 200m (0-200m, 200-400m, 400-600m... 3400-3600m, higher than 3600m), by using Model Builder in ArcGIS 10.2.

For example, if the DEM file ranges from 170-3700m, this is how I did in model builder:

I use extract by mask to clip the DEM area of interest.
I use raster calculator to select all the ranges. (Expressions: 0-200m, 200-400m, 400-600m... 3400-3600m, higher than 3600m)
extract areas of value 1.
convert files generated from step 3 to shp files.
dissolve files generated from step 4.
add a column in each shp files to indicate the class of each shp file.
assign value to each shp file.
merge all the shp files.

The attribute table of the final file would be like this:

The concept map is like this:

Is there any better way to do the classification?

I hope the model could judge the range of the input DEM file and classify it automatically, otherwise I will need to make changes every time. For instance, if another DEM file ranges from 700- 1200m. The model would have errors telling me that there is nothing exist in 0-200m and the model will stop running.


texts: Extending lines grouped by same value using ArcPy of ArcGIS Pro?

texts:
I'm trying to extend lines grouped by same value inside column.

So what i want to do is extend lines for each line with same WireID only. So for example lines with WireID 2 should all be extened to each other but not extended to lines with WireID 20.

I am trying to do it with grouping script

import itertools, arcpy

def getGrouper(row):
    """Key function to retrieve Field value. Change depending on how you group"""
    return row[0]
def getDate(row):
    """Key function to retrieve Date value. Used in mapping function below"""
    return row[1]
def mostRecentDate(group):
    """Mapping function for retrieving row in each group with most recent date"""
    return max(group, key=getDate)

# I define up here to avoid hard coding into body of script
fc = "YourFeatureClass"
# Column order is important, so define the field names outside cursor definitions
fields = ["Field","Date","Flagged_Field"]
#also putting the index of the flagged field up here, instead of hardcoding
flagindex = 2

# groups is used to store the grouped records for further operations
groups = []

# Set search cursor
with arcpy.da.SearchCursor(fc, fields) as cursor:
    # use itertools.groupby to group by Field and then append the results to groups
    for k, g in itertools.groupby(cursor, getGrouper):
        #I do nothing with the key, k, since you only care about the groups, not what they grouped on
        groups.append(list(g))

# At this point you have your groups
# The rest is to do your example of flagging the most recent date

# map groups to the most recent record in each group
flaggedrecords = map(mostRecentDate, groups)
# update records
with arcpy.da.UpdateCursor(fc, fields) as cursor:
    for row in cursor:
        # Very important to tuple the update row, or else you will not have a match
        if tuple(row) in flaggedrecords:
            # If it is in the flagged records, set the flag
            row[flagindex] = 'Y'
        else:
            # Else, unset the flag
            row[flagindex] = ''
        cursor.updateRow(row)

I came up with idea that if i calculate bearing and then add 180 to it i can use it to create extended lines. Using bearing distance to line with bearing field and new distance field (with for example 2m distance) on End_X, End_Y column will create correct extended lines. If i do the same think using 'bearing+180' field, new distance on Start_X, Start_Y column i will get other extended lines.

My logic is if i group by WireID i can start a loop that use bearing distance to line, creates intersection point and then recreate my lines using new intersection points for each WireID separately. But first i don't know how to code it and secondly is it even possible?

I am working on ArcGIS Pro Basic Licence version 1.3. Can you help me?


texts: How to use JDBC based database extension of GeoServer for other JDBC based database?

texts:
I'm doing R&D on GeoServer for my project requirement to access JDBC based Database(MapD).

So to do that GeoServer provided some extension for particular database access extension like JDBC based extension sql server,but in my scenario database is another that also have JDBC based driver.

I got github geotools link but didn't get which classes and method used for make extension and  how going on process.

So To make GeoServer extension for JDBC based database, which java classes and method need to change and make jar file to put into lib folder to access database on GeoServer.

So any help me how to do in steps.


texts: How do I export data points from topology checker into an attribute table?

texts: I used topology checker in QGIS to find all dangles in my roads polyline shapefile. The topology checker returned about 528 dangling ends, also called errors. How do I save these points into an attribute table or create a new shapefile from them?


texts:
I am having trouble getting the os.path.exists to recognize that a file exists. I keep getting "File does not exist". I reduced the code to the absolute necessary to address my question. I believe everything relevant is in there.

import os
FileName15 ="Stuck"
Folded = arcpy.GetParameterAsText(4)
if Folded == '#' or not Folded:
    Folded = "C:\Documents \test.gdb\\" # provide a default value if unspecified

if os.path.exists(Folded+FileName15+"_table.dbf") == True:
    print ("It does exist")
    print Folded+FileName15+"_table.dbf"
    # Process: Delete
    arcpy.Delete_management(Folded+Name15+"_table", "Table")

    # Process Create Table
    arcpy.CreateTable_management(Folded, FileName15+"_table",Template_Shell, "")

elif os.path.exists(Folded+FileName15+"_table.dbf") == False:
    print ("File does not exist")
    print Folded+FileName15+"_table.dbf"
    # Process Create Table
    arcpy.CreateTable_management(Folded, FileName15+"_table",Template_Shell, "")


texts:
My group and I at the University of Gent are developing a project; we made a WebApp Online, we added some layers, but now we want them to appear just when selected. How can we do this?
Also, when doing the presentation, we'd like the WebApp to appear with less zoom, I mean, show both the USA & Europe coast + UK + Iceland for example, instead of showing just the Belgian coast.
I installed ArcGIS Desktop, but still can't figure it out.
Here is the link of the webapp: http://ughent.maps.arcgis.com/apps/View/index.html?appid=3e9a6f9a2f4a42f999fc8c903076fae6


texts: I am trying to view all road names when zoomed out at this projection in QGIS - ideally at a font size of about 2. I'd then like to export this as a PDF. Can this be done?


texts: QGIS layer text field bug - unable to categorize

texts: Hi I'm still learning how to use QGIS. I've tried to add a text field called 'name' to the 'edges' layer (see screenshot attached). The layer is created but then when I create a new line on my map, the window where I would usually fill in the name seems to have been corrupted (the word 'name' has been replaced by symbols). As a result, when I try to categorize under the symbology tab, none of my categories show up and I can't differentiate the different lines I've drawn. I've tried to create the layer twice with the same issue. I'm using QGIS 3.4 long term release. Can anyone help?


texts:
I haven't setup or deployed any web services and most help seems to be around vector data... I need input on how to set up and optimize an image service with GeoServer, PostGIS to use in OpenLayers client.

Images are in 32 bit float, GeoTIFF, time-series images (Sentinel-2, Landsat-8) for disconnected areas (think like Hawaii, Montana, Ontario, Alberta), i.e. same area will have different dates of imagery over time and some areas might have images of the same date.
Want to serve old and newer imagery as they become available and so some automation will be nice. Want the user to be able to see the images and be able to view/query pixel values (e.g. NDVI).

Do I setup a folder to which images are downloaded to and this folder sync images to the raster folder in GeoServer and updates the image service? Can this be fully- or semi-automated?


texts: Setting Center Point of Pixmap Symbol for Mapserver

texts:
It seems like all the Mapserver documentation on changing the center-point of a symbol only deals with TYPE 'vector' symbols.

I have pixmap symbol (the standard looking Google marker image thingy) that i'm using with a WMS overlay. Once it's rendered the location the symbol represents ( the symbol center-point ) is the middle of the image. This makes the location the symbol represents seem off. I would rather have the marker's tip/point line up with the location it represents. I'm not sure how to do this. Any advice?


texts:
Does anyone here have experience with submitting Google indoor maps? I submitted one two weeks ago where the submission page stated "submitted for processing". (I submitted the indoor map on behalf of those I am delveloping an app for.)

I have not received any confirmation email which has worried me a little since after searching Google, I can see people in 2012 appeared to get such emails.

Is there anything I need to do beyond uploading the map? (matching the building etc. of course) - i.e. are there hidden requirements that an email or mail is sent to the address mapped or something like that? If so, I need to contact the one I am developing an app for :)


texts: Obtaining Total value of votes using spatial join in ArcMap?

texts:
I have been tasked with aggregating precinct level voting data upto the county level. I have both the layers and I used the spatial join tool to do the aggregation.

However, there is a complication. On using INTERSECT as MATCH OPTION, i dont get the total number of votes as some of the PRECINCT LINES CROSS THE COUNTY LINES.

I used several other MATCH OPTIONS but didn't get the required results.
Please suggest which MATCH OPTION should i use in order to get the total value


texts: Discrepancies in reading geopackages between ArcMap and ArcGIS Pro

texts:
I have a GeoPackage database with a few layers inside of it. One layer, which has about 1000 records, some complex geometry, but no relationships, would not open correctly in ArcMap (specifically ArcMap 10.3.1) with the error being as per screen grab (which occurs during browsing of the attribute table)

*Edit Added in the Error message as text as per Mod's suggestion.
Error Message reads "Error reading OID from table. Reading rows has been stopped.Check that the datasource is valid. Underlying DMS Error[ST_EnvIntersects: Invalid Geometry type] The operation is not support by this implementation."
However, the Geopackage is geometrically compliant inside of QGIS and other open source geometry validation tools (check validity as an example). Interestingly, ArcGIS Pro 2.2.0 has no problem reading the Geopackage.
Should I investigate the data further becuase ArcMap is telling me its invalid? Or does the QGIS OGC geometric checking tools (And ArcGIS Pro) give me enough confidence to know the geometry is valid?


texts: Using meanshift function with raster data in R

texts:
I have been playing around with different functions in R and I wanted to use the meanshift's msClustering() function with my Raster data.

Following I have a reproducible example:

a<-matrix(data=sample(1:500,25),nrow = 5, ncol = 5)
 b<-raster::raster(a)
 raster::crs(b)<-"+init=epsg:4624"
 c<-MeanShift::msClustering(a)
 b$new<-c$components

and I get an error saying

  1: In v[] <- value :
    number of items to replace is not a multiple of replacement length

which is logical given that they don't have the same length or dimensions.

So based on that example and error can anyone explain to me why I get a matrix of different dimension and give me a solution? Because I think my problem is more general and theoretical based on the classification method, rather than technical, so I figured it would be better to ask here.


texts:
I have a geopandas dataframes (A) and (B). I wanted to save them into disk using a same file but with two different layers.
I knew to save one dataframe as follows:
A.to_file('out.kml', layer='A', driver='KML')

How can I save dataframe B into same file using layer 'B'?


texts:
I am working off the OpenLayers Icon Symbolizer example (https://openlayers.org/en/latest/examples/icon.html) to create an equivalent to a Google Maps marker

new google.maps.Marker({...})

With Google Maps, I was able to add a listener and assign the ngClick event I had defined for the html element that created the marker. This way, when I clicked on the marker, it would trigger the ngClick function:

google.maps.event.addListener(scope.marker, 'click', function() {
     scope.$apply(function() {
             scope.ngClick()
     };
});

However, now that I am using the Open Layers icon, which is being added via adding a layer, I am not sure how I can link the ngClick function. This could be applied to any function that I am trying to pass on the icon click event like I had with marker. Everything I saw online, so far, has suggested that this needs to happen via map but I have not had success with it so far:

map.on('click', function(evt){
     var feature = map.forEachFeatureAtPixel(evt.pixel,
         function(feature) {
           scope.ngClick();
     });
});

While this does not throw any errors, it does not trigger my ngClick function. The example shows how to add my own popup to the icon but I want to specifically tie a predefined function to the Marker/Icon click which may be different for different markers.

Can anyone recommend what I am missing?


texts: How to process 2 GeoDataFrames by rows pairwise?

texts:
I have 2 GeoDataFrames 100 records each, and want to calculate distances between respective geometries (Points, of course).

I can't see how to make it other than write a cycle to iterate rows directly. (It seems that something changed since summer, when GeoPandas 0.4 was realeased.)

Here are the ways I used to do it:

Option 1: I assign another column. Now it sometimes turns all Points into nan.

main_df['geometry_other'] = other_df['geometry']
main_df['dist'] = main_df.apply(lambda r: r['geometry'].distance(r['geometry_other']), axis=1)
main_df.drop('geometry_other', axis=1, inplace=True)

Option 2: I join the 2 geodataframes and do the same. Same way, it turns 'geometry_other' into nans.

main_df = main_df.join(other_df, rsuffix='_other')
# then same as before

Option 3: Here's what I came up with now, after posting the question. I still don't like the repetitive .tolist() on both sides.

main_df['match_distance'] = [a.distance(b)
    for a, b in
       zip(main_df['geometry'].tolist(),
           other_df['geometry'].tolist())]

What are the options now to do this?


texts: Retire Dataset with Its Relationship Class

texts: I need to retire datasets from our SDE databases to our Retirement Geodatabase. There are several relationship classes tied to these datasets. How can I export the feature class with its relationship class intact and working. This needs to be done incase it needs to be un-retired at some point.


texts: Moving MXD to different folder without breaking links using ArcPy?

texts:
I was wondering what was the equivalent of "Save As" in ArcPy, which works perfectly fine in ArcMap. As the MXDs are currently set up with relative path moving MXD would break the links, however if MXD is saved with "Save As" option the links would be fine. I was trying to automate the procedure as there are too many files. SaveACopy (option 1) is breaking all links which is given below:

Option 1

import arcpy, os

srcdir = r"C:\temp\test1"

for basename in os.listdir(srcdir):
    if basename.endswith('.mxd'):
        mxdPathName = os.path.join(srcdir, basename)
        mxd = arcpy.mapping.MapDocument(mxdPathName)
        mxd.saveACopy(r"C:\temp\test1\test2\\" + basename)

I have tried to set the MXD to absolute path (option 2) and copy as a new copy which isn't working either (links are still broken).

Option 2

import arcpy, os

Workspace = r"C:\temp\test1"
arcpy.env.workspace = Workspace

mxdList = arcpy.ListFiles("*.mxd")

for file in mxdList:
    filePath = os.path.join(Workspace, file)
    filePathCopy = Workspace + r"test2\\" + file
    mxd = arcpy.mapping.MapDocument(filePath)
    mxd.relativePaths = False
    mxd.save()
    mxd.saveACopy(filePathCopy)
    mxdCopy = arcpy.mapping.MapDocument(filePathCopy)
    mxdCopy.save()
    mxdCopy.relativePaths = True
    mxdCopy.save()
    mxd.relativePaths = True
    mxd.save()


texts: How to define/add image id in Google Earth Engine?

texts:
This is a very basic question in Google Earth Engine.
I have an image converted from an ImageCollection using toBands():

Image = ee.ImageCollection(somelist.map(somefunction)).toBands()

Now I need to export this image to Asset. However, the metadata of this image only contains "type" and there is no "id". If I do

print(Image.id().getInfo())

the output is null. How can I specify the id by adding this property?


texts: Polygon fitting along the line keeping its area using ArcGIS Desktop (ArcPy)

texts:
I have a layer with polygons. And I need them to fit along the line. With Reshape Feature Tool I can do this, but there is one problem, the polygon area changes.
I need to keep the polygon area. Any idea how to do this?


texts: I want to mark specific lat long coordinates on a shapefile using python. I am trying to use the plot function from matplotlib to do it. In the examples that have been provided, a basemap object needs to be created first. How can i know the parameters for creating the basemap object for reading my shapefile? Also, is this the right approach to mark lat long coordinates on a map.


texts: Empty folder named '1' created when ArcMap opens

texts:
Over the past few days I have noticed a Folder named 1 on my Windows desktop (Windows 10 Home).  I delete it every time I see it, then notice it again a few hours later.

I kept a close eye on my desktop to try to catch when it appeared (I thought maybe it was a virus) and discovered that it appeared as ArcMap finished loading - the splash screen disappears and the ArcMap window appears at the same time the 1 folder appears.

No folder icon on my desktop while ArcMap is opening:

ArcMap is now open and a new folder has appeared on my desktop:

What is causing ArcMap to create this new folder on my desktop, and how do I stop it?  I have been using ArcMap for years and have never seen this before.


texts: Enabling HTTPS with GeoServer with SQL extension

texts: I have set up GeoServer (2.15.11) on a Windows server machine and have connected it to an Azure database using the SQL Server extension.  However, because SQL queries and parameters are done in plain text, I wanted to encrypt them.  Would I have to modify the SQL Server extension in GeoServer to do this, or can I implement HTTPS to do this through TLS certificates?  In the case of the latter, would I have to set that up on the server that GeoServer is on, or would that have to be done on the Azure server where the database is (which unfortunately I don't have direct access to)?


texts: AttributeError: class instance has no attribute 'class_function'

texts:
When working in my plugin's main python function, I receive an AttributeError anytime I call a class method from within the class itself.

For example the sample code below:

class PluginName:
    def __init__(self, iface):
        ...
        self.do_something()

    def do_something(self):
        ...

Would result in the error: AttributeError: PluginName instance has no attribute 'do_something'

Shouldn't do_something() be recognized as a method rather than an attribute because of the parentheses?


texts: Does Magnetic Deviation depend on altitude?

texts:
Magnetic declination, sometimes called magnetic variation, is the
  angle between magnetic north and true north. Declination is positive
  east of true north and negative when west. Magnetic declination
  changes over time and with location. As the compass points with local
  magnetic fields, declination value is needed to obtain true north.

This is the definition provided by National Center For Environmental Information.

As you can read, magnetic declination depends on location (lat/lon) and time but I would like to know if it also depends on the altitude and why it does or does not.

Note: possibly attach some reference to the answer.


texts: Calculating sum per id using UpdateCursor of ArcPy

texts:
I am trying to calculate the SUM value for (VD_SUM) Based on Basin ID and output in the field "VD_SUM_SUM".
I am not sure how to write the SUM part at the end.

This is what I have so far.

import arcpy

with arcpy.da.UpdateCursor('Streams_Geology_R1', ["BASIN_ID", "VD_SUM", "VD_SUM_SUM"]) as cursor:
    for row in cursor:
        row[2] = row[1].sum???? based on basin ID
        cursor.updateRow(row)


texts: Attribute table in QGIS showing a subset of desired features after query builder

texts:
I'm working with some OSM data in QGIS and I want to show only those points that are schools. So I go into query builder and enter "other_tags" LIKE '%"amenity"=>"school"%'.

This narrows down my features from 7000+ to only 16, which is exactly what I want.

However, when I go to view the attribute table, I only see 3 records.

I can see that there are 16 total features after the query, but 3 remain after some filter that I did not apply. I can't locate the source of this filter and don't know why it's limiting my features. In fact, I can see all the features on the map and they're all within the extent.


texts:
I am using Arcmap 10.2 and I am working with the modelbuilder.
Now I would like to create multiple maps from a single  layer file, but also use an inset map and background map. The thing is that the layer file has over 2000 records of line features. So I actually would like to cut the layer file into 100x100m (like a puzzle) and then label the line features with information based on values from the attribute table.

Can I use building map books for this and how should I do this?


texts:
I need someone who's good at GeoServer & Networks to help me figure out what I'm missing here.

I have two servers, one for database (S1) and the other one for my web app & GeoServer (S2). Both inside network A.
Application works fine on S2. Retrieving WMS and GWC layers from S1 through GeoServer [http://privateIP:8081/geoserver/wms] and [http://privateIP:8081/geoserver/gwc/service/wms]

Now, I would like to take S2 out of Network A. So I disconnected it and plugged it again to a 4G internet modem. Then changed GeoServer stores IP addresses to public IP of S1, and Layers requests to the public IP address as well: [http://publicIP:8081/geoserver/wms] and and [http://publicIP:8081/geoserver/gwc/service/wms]

At this point, requesting layers became so unbelievably slow. I can see few tiles of the layers but it takes too long to load then failed with no result but pink tiles.

I couldn't figure out what causes this. Is it rough for GeoServer to request from a remote database? and what is the suggested solution in this case?


texts: Adding new nodes in the middle of each líne of a Polygon

texts:
I am currently working on a Street Lighting Network. I have two layers in PostGIS ( PostGIS 2.1.3 and PostgreSQL 9.3.5):

Manholes (56.000 Polygon)
Pipes (Linestring)

The manholes are perfect rectangles composed by 4 nodes (vertices).

I have to digitize all the pipes from each manhole to the other manholes.

When digitizing the pipes, I have activated the snapping, but as the manholes have no middle nodes in each face I have to digitize from the vertices.

I would like to be able to draw a line between manholes but from the middle of a faces of the polygon, not from vertices.

To try to solve this problem, I thought that I could create a middle node between each of the 4 Linestrings of each Polygon. Each resulting polygon should have 8 nodes then.

My approach, which finally ended up working perfectly is:

Dump polygon into linestrings. I have found this query that does it and modify the query to construct a linestring with each startpoint, middlepoint(centroid) and endpoint. Then rebuild the polygon:

.

SELECT
id, ST_MakePolygon(ST_LineMerge(ST_Union(the_geom))) as the_geom
FROM(
SELECT id, ST_MakeLine(ARRAY[sp,st_centroid(ST_MakeLine(sp,ep)),ep]) as the_geom
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT id,
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
       -- extract the individual linestrings
      (SELECT id, (ST_Dump(ST_Boundary(the_geom))).geom
       FROM manholes
       ) AS linestrings
    ) AS segments
) as lines
GROUP BY
  id

Result:


texts:
Is there any easy way to resample 30 meter NDVI from Landsat 8 images to 5 metres?

I use raster to point and the interpolation of point to raster using Kriging method.

I know that it may not be a proper way of  getting 5 m NDVI but for the purpose of demonstration I need to do.


texts:
I'd like to add the geopdf format as output of a library that I am working on. After have looking around, I have not been able to find a (c, c++, python) library that meet my license requirements  (LGPL, BSD or MIT). The only way that I found until now is GDAL compiled against libpodofo as reported here. However, it looks a little bit cumbersome:

  However, for getting the imagery, the pdftoppm utility that comes with
  the poppler distribution must be available in the system PATH.

Do you know any alternative solution?


texts: Upgrade Geoserver version running with Jetty

texts:
I have a Geoserver install (2.1.3) running with Jetty that I need to update to the latest release, but I can't lose any of the residing data files and layers configuration. Any help to obtain this (on an Ubuntu 10.04LTS VM).
Thanks.


texts:
I have a table of point geometries like pointID, village and geom. I want to group the point by village and create a centeroid representing the village.

How do I write that query in PostgreSQL?

I am using PostgreSQL 9.6 with PostGIS 2.2.


texts:
I have created a contour in QGIS based on a TIN interpolation.
The contour lines are displayed over the whole area of the TIN. Is there a way to limit the display to the layer of the outer border of the lot?


texts: Sample all existing points of a polygon in Python

texts:
I'm trying to sample raster values to points based on the geometries of a shapefile using the below code but the raster values are only sampled at the coordinates of the polygon geometry. Is there a way to sample for all the points along the boundary and inside the polygon.
multi_values_points = pd.Series()
 with fiona.open(inputShape) as shp:
     for feature in shp:
        Name = feature['properties']['Name']
        coords = feature['geometry']['coordinates']
        with rio.open(inputRaster) as stack_src:
           value = [v for v in stack_src.sample([coords])]
        multi_values_points.loc[Name] = value


texts: QGIS: How can I reproject multiple associated layers when improving the georeference of an ancient map?

texts:
If I add or move some georeferencing control points and re-georeference a map, then is there a way in which I can make the latest QGIS reproject the points/lines/polygons in associated layers automatically?

I georeferenced a 17th-century map, then digitized thousands of map elements (points, lines, polygons) in various additional layers.  Now I find myself wanting to add or move some georeferencing points to improve the location of features in parts of the map, which means that the corresponding digitized map elements (points, line vertexes, polygon vertexes) in other layers have to move, too.  I did this once by hand: re-georeference the map, then manually move all points/lines/polygons as needed.  This was a lot of work.  How can I have QGIS do most or all of that work?


texts: Google Earth land features importing to ArcMAP

texts:
I'm wanting to import land features, such as lakes, reservoirs and other land-made features onto my ArcMAP. Is that possible as this is new with my client wanting this done.

What needs to be done as I've been researching this.


texts: How to evaluate the hilliness of an area using SRTM data?

texts:
I want to evaluate the hilliness of an area (e.g. municipality) using SRTM data. I thought about just calculating the standard deviation of the raster cells but my gut feeling is that this approach might be to simplistic.
What is the best approach to measure/evaluate the hilliness of a region using raster data?


texts: plm/splm in R - fixed effects panel models with time invariant variables

texts:
a bit new to working with panel data models. i'm estimating the following model for crime rate:
tot_crime <- log_rt_crime ~ ndvi + log_disad_inde + log_density
where ndvi is a time-variant continuous variable, and log_disad_inde and log_density are time-invariant continuous variables. in brief, here is the data structure:
geoid           log_density     log_disad_inde  time    ndvi    log_rt_crime
170310101001    0.6633400        0.9688232      1   0.21318862   2.213387
170310101001    0.6633400        0.9688232      2   0.46648943   1.961499
170310101001    0.6633400        0.9688232      3   0.39554212   2.725013
170310101001    0.6633400        0.9688232      4   0.14324121   2.414422
170310101002    0.8935242        1.1857459      1   0.13986452   3.758296
170310101002    0.8935242        1.1857459      2   0.26278001   3.559268
170310101002    0.8935242        1.1857459      3   0.20211675   3.646066
170310101002    0.8935242        1.1857459      4   0.09324110   3.390321

when i estimate the individual fixed effects within estimator without spatial effects:
summary(FEwithin_nosp_crime <- plm(tot_crime, data = final_long_panel, model = "within", effect = "individual")),
R understandably only estimates my time-variant variable, ndvi, and drops log_density and log_disad_inde. however, when i account for spatial dependency,
summary(FEwithin_SAR_crime <- spml(tot_crime, data = final_long_panel, listw = queen, lag = TRUE, model = "within", effect = "individual", spatial.error="none")),
i do get estimates for my time-invariant variables. why is this the case? cheers.
Note: above model is the spatial lag model; i also am getting estimates for time-invariant variables when i estimate the spatial error model and other spatial panel specifications.


texts: Looping through list of coordinates to create polygons on Google Earth Engine

texts:
I have imported a CSV file with just one column. Each row of the column has multiple coordinates required to draw a polygon. (This was collected on the ground)

I have used the below code to extract the coordinates from the CSV as a list

var poly = table2.aggregate_array('wkt')  //wkt being the column name

Now essentially, I would want to loop through the list, create and display the polygons on the map.

I am looking for something like this:

var newPolygon = ee.Geometry.Polygon(poly[i])

If there is a simpler way to use the coordinates to plot polygons can you let me know?


texts: Arcobjects in Python, inherit from Abstract Class

texts:
As for question, I am trying to adjust Grid Map labels via Arcobjects with Python. For now I do have:

from comtypes.client import CreateObject, GetModule
import arcpy

def CType(obj, interface):
   """Casts obj to interface and returns comtypes POINTER or None"""
   try:
       newobj = obj.QueryInterface(interface)
       return newobj
   except:
       return None

def NewObj(MyClass, MyInterface):
   """Creates a new comtypes POINTER object where\n\
   MyClass is the class to be instantiated,\n\
   MyInterface is the interface to be assigned"""
   from comtypes.client import CreateObject
   try:
       ptr = CreateObject(MyClass, interface=MyInterface)
       return ptr
   except:
       return None

esriCarto = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.0\com\esriCarto.olb")
esriCartoUI = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.0\com\esriCartoUI.olb")
esriMapUI = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.0\com\esriArcMapUI.olb")
esriFrame = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.0\com\esriFramework.olb")

arcpy.SetProduct('Arcinfo')

pApp = NewObj(esriFrame.AppROT, esriFrame.IAppROT).Item(0)
pDoc = pApp.Document
pMxDoc = CType(pDoc, esriMapUI.IMxDocument)
pLayout = pMxDoc.PageLayout
pGraphContLayout = CType(pLayout, esriCarto.IGraphicsContainer)
iFrameElement = pGraphContLayout.FindFrame(pMxDoc.ActiveView.FocusMap)

As far as I understand, iFrameElement is an interface of an abstract class from which I need to inherit attributes (pointer) to MapFrame object. How do I do that? How do it get to object with IMapGrids interface. Any suggestions?


texts: Show wms feature infos after map is rendered (OpenLayers)

texts:
I display wms features on a map and i want to see the attributes of all features after the map is loaded.
I know the GetFeatureInfo request and i made examples when i clicked on a feature and than a popup was appeared with the infos. But i need a solution when all the popups appear after the map and the features are loaded.
Any idea?


texts: Connect ArcMap to PostgreSQL database on other than default port

texts:
I have a PostgreSQL database on a computer, hosted on port 5433. on the same computer I also have a database hosted on port 5432

I am trying to connect to the database via ArcMap database connection. typically this is pretty straighforward (if the pg_hba.conf is setup to allow connections)

simply enter the computer name in instance and put in the sign in and password credentials.

I do not know how to connect to the database on port 5433 from another computer


texts: Create groups according to spatial contiguity?

texts:
I have a large set of points for school sites and I want to divide these schools into 6 groups, for example, and I want each group to consist of 200 schools, provided they are spatially close to each other.
How can I do that in ArcGIS Pro?


texts: Combine bands from two image collections

texts:
I am making a image collection of NDVI from Landsat 5 over a particular location.  I want to run a regression on each pixel through the NDVI using the per pixel NDVI value as the dependent variable and the year of image acquisition as the independent variable.

So far I am doing this:

// Define start and end dates and geometry to select images by
var start = '1984-05-01';
var end = '2011-09-30';
var polygon = ee.Geometry.Polygon([[
 [-97.49404907226562, 46.59473135600069],
 [-97.09442138671875, 46.59567501063883],
 [-97.05734252929688, 46.31184150036163],
 [-97.53524780273438, 46.30709840788667]
]]);

// function which will add a band to original image stack based on start year
function createTimeBand(img) {
  var year = ee.Date(img.get('system:time_start')).get('year').subtract(1984);
  return ee.Image(year).byte().addBands(img);
}

//Get an image collection of images of interest and add the band from function created above
var images = ee.ImageCollection('LANDSAT/LT5_L1T_TOA_FMASK')
.filterDate(start, end)
.filter(ee.Filter.dayOfYear(120, 275))
.filter(ee.Filter.eq('WRS_PATH', 30))
.filter(ee.Filter.eq('WRS_ROW', 28))
.filter(ee.Filter.lessThanOrEquals('CLOUD_COVER', 10))
.filterBounds(polygon).map(createTimeBand);

//make ndvi with a function
var make_ndvi = function(image) {
   return image.normalizedDifference(['B4', 'B3']);
 };

var ndvi = images.map(make_ndvi);

//get an image collection of the bands added which represent years since start of time series

function get_constant(images){
  return images.select('constant');
}
var constants = images.map(get_constant)

This code results in two image collections with 103 images each. The collections are namedndvi, and constants and each have 1 band per image,  I want to combine the individual bands from the two collections in a sort of zipped list (in python its called a zip list at least) so that my output is one image collection with 103 images and with two bands per image.  Then I can use regression on the bands through the collection. This problem is similar to the linear_fit code in the examples of Earth Engine, but I have to modify their code since I lose the original metadata when making NDVI.  If I can keep the metadata when making my NDVI then I can solve it from there too.


texts: Which data formats are best for Open Data distribution?

texts:
What are the pros and cons of different data formats (performance, file size, etc.) when considering open data distribution?

Our organisation wants to publish data as open data. However, there is no clear idea on which data formats to use. Ofcourse, the more 'open' a data format is, the easier it is to use.

Which data formats are the most 'open' and therefore most usable for the distribution of Open Data when taking the following types in consideration?:

raster data (I'm thinking: GeoTIFF, Erdas Imagine IMG?)
vector data (I'm thinking: GML, CSV, ESRI Shapefile, DXF?)
tabular data (I'm thinking: CSV?)
3D data (I'm thinking: CityGML?)
3D point coulds / LIDAR (I'm thinking: LAS?)
am I forgetting something here?

Also, if there is documentation about open data formats I'm very interested if you would like to share.


texts: Remove holes from polygons using geometry objects

texts:
I've written a script that should remove holes from polygons with use of cursors and geometry objects. However, holes remain even though I'm constructing new polygons from the vertices from the first part of the input polygon geometry.

Here's the script:

def RemoveHoles (inFc, workspace):
    outFc = UniqueFileName(workspace) #creates a file path/name string
    array = arcpy.Array ()
    sr = arcpy.Describe (inFc).spatialReference
    outPath, outName = os.path.split (outFc)
    arcpy.CreateFeatureclass_management (outPath, outName, "POLYGON", spatial_reference = sr)
    with arcpy.da.InsertCursor (outFc, "SHAPE@") as iCurs:
        with arcpy.da.SearchCursor (inFc, "SHAPE@") as sCurs:
            for geom, in sCurs:
                for i in range (geom.partCount):
                    part = geom.getPart (i)
                    for pnt in part:
                        print pnt
                        if not pnt:
                            break
                        array.add (geom.getPart (0))
                    print array.count
                    polygon = arcpy.Polygon (array)
                    array.removeAll ()
                    row = (polygon,)
                    iCurs.insertRow (row)
    del iCurs
    del sCurs
    return outFc

For this polygon:

The output reads as so:

571048.0859 4230953.3243 NaN NaN
571762.4623 4225079.5626 NaN NaN
566206.2012 4227698.9428 NaN NaN
571048.0859 4230953.3243 NaN NaN
None
4

And yet, despite the fact that I'm constructing a new polygon from three points, the hole is still preserved in the output feature class. Does anyone have any idea why?


texts: Extracting vector/point data from single-layer, non-georeferenced vector PDF file

texts:
I am trying to extract point data from a non-geo-referenced vector PDF file, to import into QGIS/ArcGIS. You can download a copy of the PDF file here, if interested, but I've tried to include screenshots of the relevant parts to make that unnecessary since it's such a large file (~171 MB).

The PDF file is a map of mineral resources in Afghanistan, which looks like this:

Here is a closeup image of a small segment of the map -- notice that it contains a bunch of labeled locations (cities/provinces, etc), colored regions (representing geological/rock types), and point data with different colored/shaped icons (representing known mineral/petroleum deposits):

I am primarily interested in extracting the locations of the little icons representing the geo-locations of mineral/petroleum deposits, and using them to generate point data that I can import into QGIS. That is, for all of the icons listed in the legend below, I want to find any instances on the map and create a table that has their x,y location and "category" field that matches the categories listed in the legend:

The PDF file does not appear to be geo-referenced, but the map is drawn to scale and surrounded by a lat/long-grid. The map contains the following information about the scale (and authors):

In the bottom left corner there is the following information about the projection which I feel could somehow be used to map from (x,y) coordinates in the PDF file to (x,y) coordinates in WGS84:

Obviously, if there is an easy way for me to do this in QGIS, I'd love to know about it ... but from my searches so far, it's appears likely that I'm going to have to write code to extract the data from the PDF ...

If that is the case, then I would be interested in Python/C/C++/Scheme libraries that could be used to extract the information I described above. And if no libraries exist that can do this easily, then please share ideas about resources/techniques that might help me code my own solution for this.

Any ideas on how to tackle this problem would be useful, even if they only solve some small part of it ...

UPDATE:

So I've experimented a bit with a FOSS tool called pstoedit which can extract some of the data from the PDF file and convert it into a DXF file (a common CAD format that QGIS can work with). I used the DXF with splines format option (pstoedit -f dxfs infile.pdf outfile.dxf) to convert the file, and got a DXF file that QGIS was able to import using WSG84 (I still haven't figured out how to use the projection/reference info in the bottom left of the map to overlay this point data, in the right locations, over my base map of Afghanistan).

However, there is very little organization to this data - it's just (x,y) data and generic Text field names like "Text" or "EntityHandle" (see left sidebar in image below). And as user30184 pointed out in their answer, all of the data is jumbled together in one Layer (0). So there is no clear way for me to associate the graphical icons with the individual text labels (mineral deposit categories) from the legend. Here is what it looks like so far - I'll post more if I make more progress:


texts: extract() function from r raster package returns different cells if "weights" parameter is TRUE or FALSE

texts:
I am trying to extract cell values from a RasterLayer based on a SpatialPolygons object. I create the raster with

library(raster)

rasterValues <- matrix(rnorm(20),4,5)
r <- raster(rasterValues)
extent(r) <- extent(13.875,14.5,45.125,45.625)
projection(r) <- crs("+proj=longlat +datum=WGS84")

and then I create the SpatialPolygons with

x <- c(14.00, 14.27, 14.27, 14.00)
y <- c(45.25, 45.25, 45.5, 45.5)
spPoints <- cbind(x,y)
sPoly <- Polygon(spPoints)
sPolys <- Polygons(list(sPoly),1)
sPolyg <- SpatialPolygons(list(sPolys), proj4string = crs("+proj=longlat +datum=WGS84"))

As you can see from the picture, the polygon (colored in black) totally covers 4 cells and a small part of two cells of my raster.

When I try to extract values with weights=TRUE, I get:

cellValue <- extract(r,sPolyg, cellnumbers=TRUE, weights=TRUE)
cellValue <- cellValue[[1]]
cellValue

      cell      value     weight
[1,]    7 -0.1658694 0.22727273
[2,]    8 -0.8996764 0.22727273
[3,]    9 -0.6951724 0.04545455
[4,]   12  0.5717841 0.22727273
[5,]   13  1.5263394 0.22727273
[6,]   14 -1.2750113 0.04545455

When I try to run the previous command with weights=FALSE, I get:

cellValue <- extract(r,sPolyg, cellnumbers=TRUE, weights=FALSE)
cellValue <- cellValue[[1]]
cellValue

cell      value
[1,]    7 -0.1658694
[2,]    8 -0.8996764
[3,]   12  0.5717841
[4,]   13  1.5263394

The two cells that are only partly covered by the SpatialPolygons object are not returned.

The question here is: why are some cells not returned, when we set the "weight" parameter to FALSE?


texts: How to publish a service on ArcGIS Server using versioned data in ArcGIS 10.3?

texts:
Fairly new to GIS here, especially new to working with enterprise geodatabases. I am attempting to publish a service from arcmap to arcgis server using a version of the DEFAULT. I am using connection files from two separate databases in SQLserver 2008 (db1 and db2). Let me walk you through my processes:

add database connection using sde database authentication to db1
import data, manage privileges, register as versioned, add Global IDs
geodatabase administration > create new version from DEFAULT (named v1)
add new database connection using my windows authentication to db1
change geodatabase connection properties > select v1 from Transactional version
add data to arcmap from v1 version

Register database
When I register my geodatabase with ArcGIS Server, I import the connection file to db1 as the publisher database connection, and the connection file to db2 as the Server database connection, leaving the create geodata service option unchecked. The database registers and validates properly on arcgis server manager. The reason I have the separate server database connection is so I can replicate my data and synchronize my edits between my v1 version and arcgis online.

Next:
- right click my windows authenticated geodatabase > share as geodata service (so i can publish it in a specific folder)
- go through the publishing steps but when I click analyze, I keep getting the error 00133 Geodatabase is not registered with the server.

HOWEVER, when I follow these same steps using the DEFAULT version, it publishes fine. Arc does not seem to produce a separate connection file for the versions I create, so how do I get my service to recognize my registered geodatabase when publishing a version???

The same thing happens when I publish a map service (with feature access enabled) using v1 versus DEFAULT, except I get the error 00090 Feature service requires a registered database.