texts: How do I determine the center point of several features in OpenLayers?

texts:
I have several polygons being drawn on a map, usually representing a country or province or state.  Depending on some chosen data, different regions may be selected, so what I was hoping is that there's a way to determine the center point between all of the features on the layer so the map can automatically be focused on the most relevant point.
any help would be great.
Thanks


texts: Draw polygon with gradient color in Geoserver

texts:
I have polygon like sector. I want to fill it with gradient color.

Is it possible to do this in Geoserver?

Thanks.

UPDATE

I found some code with filling image

<FeatureTypeStyle>
 <Rule>
   <PolygonSymbolizer>
     <Fill>
       <GraphicFill>
         <Graphic>
           <ExternalGraphic>
             <OnlineResource
               xlink:type="simple"
               xlink:href="colorblocks.png" />
             <Format>image/png</Format>
           </ExternalGraphic>
         <Size>93</Size>
         </Graphic>
       </GraphicFill>
     </Fill>
   </PolygonSymbolizer>
 </Rule>

i can put here big circle which is gradient, but does it has rotation or position of image?


texts: Getting household income data based on address to block level?

texts:
I am beginning a research project and need to find a way to access the census data about household income.  We need it at a closer level than by zip code since income varies a lot by zip code in our area.

Is there a way to get household income based on address or at least down to the block level?


texts: Is it possible to make a routable graph from Polygons?

texts:
I would like to know if it is possible to make a graph with pgrouting 1.05 on polygons. What I mean is that I have to make a graph on a database without roads but with polygon geometry. I would like to make a graph on this geometry.

Thanks for your help and I'm sorry for my english


texts: Accessing a non spatial table in an ArcGIS Map Service in QGIS

texts: I have a client who wants some information that is stored in a non-spatial table on an ArcGIS Map Service.  When connecting to that service via either the Map Service Connection or Feature Service Connection in QGIS the table does not show up.  Is there a way to get at that table in QGIS without knowing the architecture of their server or do I need to have my client ask the city to make the table part of their open data website?  It's the table located here https://gis.burnaby.ca/arcgis/rest/services/BurnabyMap/BBY_PUBLIC_TBL/MapServer/1


texts: Adding line from one long/lat point to another long/lat point using CartoDB

texts:
I have imported an excel s/sheet into CartoDB.

I have the following columns:

BookingReference  | Long | Lat

I want to see the journey my app users take. As such I want add a line between long lat points for a booking reference.

For example:

0001/34458      51.5074 -0.1278
0001/34458      52.8984 -1.2698
0001/34458      53.8008 -1.5491
0001/34459      21.1619 -86.8515
0002/44492      21.1619 -86.8515
0002/44492      51.5074 -0.1278

One booking reference might have several long/lat points, others may just have a single one or two+


texts:
I'm planning a research that includes couple Ikonos, QuickBird and GeoEye digital images.

Until now, I have had no luck getting answer from the official website of Landinfo.

Is there any other way to order this imagery, and is there a minimum order price (e.g. minimum $300)?

I also cannot see any 'Buy Now' options on their website.

Are the orders strictly via e-mail?


texts: How to Spatial Join centroids of parcels to point data and get the distance

texts:
I'm currently using ArcMap 10.1 (basic license). I have a LSOA shape file of Camden borough. From there I added 2 extra fields within the attribute table, X and Y. I then calculated the 2 extra fields to give me the centroid of each output area. I also have a point data shape file of the underground stations within Camden borough. three things I need to work out:

Get the distance for every LSOA to the closest underground stations in order to create a choropleth map
I need to create a polygon vector file containing the catchment areas of every underground station in Camden. If there are 19 stations then I assuming that there are 19
 catchments areas
Get the centroid of the catchment areas and then get the distance to the closest underground stations and create a choropleth map

I hope that makes sense


texts: Pixel value in export GeoTIFF file stretched by color value

texts:
I have exported GeoTIFF files from GEE and attempted to keep the pixel value in that raster layer. Once the file was imported into ArcGIS Pro the stats (min, max, mean and sd) were stretched between 0-255 instead of keeping the original pixel value from NDVI.
Below is my code:

//get mean and std of NDVI
var mean = ee.Number(stats.get('NDVI_mean'));
var std = ee.Number(stats.get('NDVI_stdDev'));
//find regions in ndvi has value higher than mean +1SD
var maskImage = ndvi.updateMask(ndvi.gt(mean.add(std)));
//use the visualize function to set visualize parameters min and max equals to min and max in the raster
var visualizeMinMaxNDVI = function(maskImage) {
  var minMax = maskImage.reduceRegion({
  reducer: ee.Reducer.minMax(),
  bestEffort: true,
});
  var visParams = {
          min: minMax.getNumber('NDVI_min'),
          max: minMax.getNumber('NDVI_max')
          };
  return maskImage.visualize(visParams)
            .set({min: minMax.getNumber('NDVI_min'), max: minMax.getNumber('NDVI_max')});
};

// set visualization parameters based on min and max value
var maskImage = visualizeMinMaxNDVI(maskImage);

// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: maskImage,
  description: '21072016UNDVItest',
  scale: 10,
  region: table,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});

I'm very new to GEE, I think when I export my GeoTIFF layer the data was 8 bit (shown in ArcGIS Pro) which caused the color stretch?


texts: Most efficient way to retrieve pixel distributions inside many polygones (custom zonal stats)

texts:
I want to get the distribution of pixels from a 1 band raster that fall in each polygon of a shapefile to proceed to some statistics. I guess a zonal statistics does it in some kind of way but I need other statitsics that the ones provided by the tools.

Specifically, I want to sort the pixels inside the polygons and retrieve the 50th value (if exists).

Converting my raster to points and joining them spatially to my polygons and exporting my table for the statistical analysis would do it but I will obtain more than 80M points and I estimate that it will take more than 4 days if no issues.

Does anyone know how to perform this in a more efficient way ?


texts:
I have run the gdaltindex command with output to a MySQL table, however only a list of file name are

gdaltindex -lyr_name Source_raster_idx -f MYSQL MySQL:mapDB,host=localhost,user=admin,password=pass */*/*/*/*.jp2

It loads the filename and path into the column named "location", However, I don't see where the spatial coordinates are supposed to be stored.

Am I supposed to add a column for this?


texts: How to display ALL raster output from batch processing in ArcGIS 10.4

texts:
I am batch converting ASCII files into ArcGIS geodatabase using a simple ASCII to raster model. I was using ArcGIS 10.3.1 and everything was running great, all my surfaces/files were showing up at the end of the batch processing in my table of content/display with the proper symbology etc.

My company upgraded to ArcGIS 10.4 and now only the last raster in my batch list shows up in my display despite having all the options turned on for "Add result to display" (in both my model and in Geoprocessing Options).

Here is a screenshot of my model, output raster parameters, Geoprocessing parameters, batch window and resulting map window:

I only have two files in this exemple but I'm usually working with with tens of new files every day.

My input file is in ESRI ASCII format.

I did all possible combinations of "Add to Display" options i.e. only the first output, only the last output, all three.

But in the end only the last file in the batch list gets displayed.


texts: Merging duplicate features in a dataset but keeping all the attributes

texts:
I have a polygon shapefile, in which there are many duplicate records (duplicate in the spatial sense, but with different attributes). What I want to do is merge all features that have the exact same geometry but keep the attributes from all of the polygons. Is this possible to do?

For an example, if I had the following as an input and all three features has the same geometry.

Then I want something like this at the end


texts:
I'm currently working on a project based around cycling. I need to know which dataset I should get hold of which will allow me to create a multimodal network in Network Analyst so that when I run the network, cycle paths such as the National Cycle Network are also utilised by bikes. I have been using OSM data for networking, but i'm unsure if this can be integrated with a cycle network.

What dataset would provide me with the information I require?


texts: I've been trying to find out what the projected bounds are for EPSG:3575 but have not been able to find a definitive answer. Normally, I'd just use the values from EPSG.io but when I tried these in mapserver and geowebcache they both threw back error messages saying 'min x, miny, max x, max y not valid'. I've also tried the values from spatialreference.org and whilst these are valid, the cache gridset I created in GWC then did not work with several WMS services from other organisations. I have also found this site http://nsidc.org/data/atlas/epsg_3575.html, but this just gives the same info as the EPSG registry in effect. Ive tried converting the 4326 bbox to 3575 in a couple of programs (ogr and a web conversion site), but these just out with differing answers, neither of which look correct. any help would be gratfully received!


texts: Getting maximum extent of 2 different layers for map display using PyQGIS?

texts:
In the QGIS 3 layout, I would like to display two layers and for this I need their combined maximum extent. So far I've got:

import os
import processing
from qgis.core import *
from qgis.gui import *
import qgis.utils
import shutil
from qgis.PyQt import QtGui

project = QgsProject.instance()
project.setCrs(QgsCoordinateReferenceSystem(21781))

iface.mainWindow().blockSignals(True)

# add vector layer
shplayer = iface.addVectorLayer("C:\\MyData\\A300.shp", "A300", "ogr")
shplayer.setCrs(QgsCoordinateReferenceSystem(2056))
project.addMapLayer(shplayer)

# add raster layer 1
rlayer1 = iface.addRasterLayer("C:\\MyData\\Area_1.asc", "Area1")
rlayer1.setCrs(QgsCoordinateReferenceSystem(2056))
project.addMapLayer(rlayer1)

# find max extent of all layers
ext_r = rlayer1.extent().toString().replace(":",",").split(',') #(xmin, xmax, ymin, ymax)
ext_shp = shplayer.extent().toString().replace(":",",").split(',')
ext_all = [min(ext_r1[0],ext_shp [0]), max(ext_r1[1],ext_shp [1]), min(ext_r1[2],ext_shp [2]), max(ext_r1[3],ext_shp [3])]

#-- Create custom caneva
manager = project.layoutManager()
layout_name = 'my_custom_layout'
layout_list = manager.printLayouts()

layout = QgsPrintLayout(project)
layout.setName(layout_name)
layout.initializeDefaults()
manager.addLayout(layout)

#--- create map item (shapefile, raster 1)

map1 = QgsLayoutItemMap(layout) # create Layout item
map1.setRect(QRectF(10, 10, 210, 90)) # set the map item rectangle
map1.setPos(10,10) # map item upper left corner position
map1.setFrameEnabled(True) # to have a map black border

ms = QgsMapSettings()
ms.setLayers([shplayer, rlayer1]) # set layers to be mapped

# !! here I'd like the max extent of both layers! QgsRectangle(rlayer1.extent() + shplayer.extent())
map1_extent = QgsRectangle(float(ext_all[0]), float(ext_all[2]), float(ext_all[1]), float(ext_all[3])) # must be (xmin, ymin, xmax, ymax) but does not work!
map1_extent.scale(1.0)
ms.setExtent(map1_extent)
map1.setExtent(map1_extent)

layout.addLayoutItem(map1)
map1.attemptMove(QgsLayoutPoint(10,10))
map1.attemptResize(QgsLayoutSize(210,90))

# export map in .pdf
exporter = QgsLayoutExporter(layout)
exporter.exportToPdf("C:\\mapstest.pdf", QgsLayoutExporter.PdfExportSettings())

project.write()

How can I fix this?


texts: Setting Reactivity in Leaflet Shiny to subset and redraw features from two fields

texts:
I am new to Leaflet and Shiny and would like to create two reactive functions from two columns in my shapefile frame(map1). I would want the application to redraw polygons based on the input made in these functions (input of an Area and an Integer value).
How do I fix the below code to derive the desired output?
ui <- fluidPage(
    titlePanel(p("Learn Shiny", style = "color:#3474A7")),
    sidebarLayout(
        sidebarPanel(

            checkboxGroupInput(inputId = "AreaFinder",
                               label = "Select Areas(s):",
                               choices = unique(map1$Areas),
                               selected = "Christchurch"),

            checkboxGroupInput(inputId = "NumberID",
                               label = "Select number:",
                               choices = unique(map1$No),
                               selected = 0)
            ),

        mainPanel(
            leafletOutput(outputId = "Dorset_map")
        )
    )
)

# server()

# server()

server <- function(input, output) {

    data_input <- reactive({
        map1 %>%
            filter (map1$No == input$AreaFinder)
        AND
        filter(map1$Areas == input$NumberID)
    })

    output$Dorset_map <- renderLeaflet({
        leaflet() %>%
            addTiles() %>%
            addPolygons(data = map1,
                        fillColor = pal(data_input()$Area),
                        weight = 1,
                        smoothFactor = 0.5,
                        color = "white",
                        fillOpacity = 0.8,
                        highlight = highlightOptions(
                            weight = 5,
                            color = "#666",
                            dashArray = "",
                            fillOpacity = 0.7,
                            bringToFront = TRUE))

    })
}

# shinyApp()
shinyApp(ui = ui, server = server)


texts: QGIS - How to add joined values table data to shapefile

texts:
For some reason I lost my previous question here but I proceeded a little bit after upgrading to the newer version of QGIS.

My problem was combining spreadsheet .csv data to the .shp map layer, where the joining and target field were building values and it only recognized the first feature and missed all the other features coming from the same building. I managed to do a join with join attributes by field value (one-to-many) and the joined table now includes all the features.

Now I'm trying to join this new joined table to building point data map layer .shp, but it does the same thing again leaving the duplicate values out (building value). How can I add this joined table with all the features in it to the point data map layer .shp? Should I use different join values? I can only join the spreadsheet and building data with a certain building value.


texts: ArcObjects - IGraphicsContainerEvents_ElementAddedEventHandler Unregisters After a Few Fires

texts:
I am using IGraphicsContainerEvents_ElementAddedEventHandler to do stuff with graphics added by users in an ArcMap session using the standard drawing toolbar.  As the code snippet below shows, its a pretty straight-forward process, and works perfect for the first few (2-3) graphics added.  However, eventually 'ElementAdded' stops firing for no apparent reason (i.e no errors that I can catch).

Does anyone have any idea of what would be unregistering the listener (OnElementAdded) and/or preventing the ElementAdded event to be fired?

ICompositeGraphicsLayer activeGraphicsLayer = (ICompositeGraphicsLayer)mxDoc.ActiveView.FocusMap.ActiveGraphicsLayer;
IGraphicsContainerEvents_Event containerEvents = (IGraphicsContainerEvents_Event)activeGraphicsLayer;
containerEvents.ElementAdded += new IGraphicsContainerEvents_ElementAddedEventHandler(OnElementAdded);

private void OnElementAdded(IElement Element)
{
// do stuff
}


texts:
I'm looking to make a Google Earth Engine App that automatically updates on a regular basis as new satellite imagery comes in for a given dataset. There doesn't seem to be a way to do this easily from within the GEE code editor / 'publish app' menu – is it possible to interact with an already-published GEE App via the GEE API?
One typically specifies a time frame for resources used in the GEE code editor, e.g.:
var dataset = ee.ImageCollection('TRMM/3B42')
                  .filter(ee.Filter.date('2019-03-15', '2019-03-16'));

For automatic updates, the updated time & date would need to be specified in the image query code.
Automatic updates would be useful to create, for instance, an app that displays latest weather changes / cloud cover movement.


texts: Is there a way to save the georeferencing residuals table in ArcGIS Pro?

texts: I am conducting a georeferencing task for which I need to calculate statistics on the x and y residuals. I note that these are helpfully calculated in the control point table both before and after georeferencing Is there a way to copy or save these values from these tables? I note that it is possible to save the original control points, without residuals, but saving after the georeferencing does not seem to change the coordinates from the originals.


texts:
I have a Vector layer containing 36000+ line Objects with a large number of attributes that is currently saved as Shapefile on a server (local network).

Multiple users access and edit this Layer from different computers in the local network via the attribute table. The workflow seems to get slower  as soon as more than one person works with the layer.

Is a shapefile the right solution for multiuser access?
Where to save the Layer and in what format?


texts: Creating new polygon in between other polygons?

texts:
I want to fill the missing space and create a polygon out of it. How do you do that with QGIS?
The missing space is the white part between the selected/yellow polygons.
Here is a snapshot:


texts: How can I fix QGIS 1.8.0 Spit Plugin "connection failed"?

texts:
I want to import shp to to postgreSQL with the spit plugin. Testing the database connection to my PostGIS database ist successfully (I already imported layers with this connection via the "add postGIS Layer" button), but if i try to connect to my db within the spit dialogue, entering the password, there will be an error message "connection failed".

I tried to use the DB Manager, three times I entered the password, then error message "unable to connect fe_sendauth: no password supplied" .. don't know, if this error is related to the above.

Have got Win7 with QGIS 1.8 and PostgreSQL 9.2. I supplied a passowrd during PostgreSQL installation to the superuser.
Editing pg-hba.conf without success. In fact, couldn't connect in pgAdmin after appending the conf-file.
hopefully, going to do the whole QGIS-PostgreSQL-installation on another machine this month.

Thanks to your answer and comment, @underdark und @RomaH!


texts: Encoding error when importing to PostgreSQL using osm2po SQL file

texts:
After running osm2po on a pbf file for Europe I am importing the resulting files into a PostGIS and pgRouting enabled database. During the import psql it reporting many encoding errors because the WIN1252 found encoding has no equivalent in UTF8.

How can I solve that issue?

The initial pbf file was downloaded from GeoFabrik and osm2po used in the latest version available on Anitas webpage. PostgreSQL is in version 9.2, pgRouting is version 2.0


texts: I need to access tables which are in an Oracle database by ArcCatalog. I am new to this ArcMap software, so please suggest to me detailed steps to access database tables from ArcMap.


texts:
I have been using PostGIS 3.0. The ST_HausdorffDistance() function calculates a symmetric distance between two polygons. However, Hausdorff distance is not a metric since it is the maximum distance between the closest pairs of vertices. Therefore, it depends on the number of vertices each polygon possess. I could not understand how PostGIS made it symmetric.
In addition, when I investigated a toy dataset, where we have two polygon layers (A and B), and their corresponding vertices generated by the Coord. AttribuTable plugin of QGIS (A_vertices, B_vertices), I calculated the Hausdorff distance from A to B as 16.41 and 18.13 for polygons with IDs 111 and 222 respectively. These distances are 31.19 and 25.70 in PostGIS.
select ST_HausdorffDistance(a.geom, b.geom)
from a, b
where a.id=222 and b.id=222


texts: EasyCustomLabel in QGIS do not resize the labels

texts:
I use EasyCustomLabel to move some labels, but if I try to change the dimensions of the chars I cannot view any changes, both in the work and in the printing sheets.

I use QGIS 2.6.1.


texts: st_3dshortestline equivalent in Python 3

texts:
I have two geopandas Dataframes with 3D geometry. I could transform it in any other data type it could be more easy to use.
I would like to know if there is something similar to the st_3dshortestline PostGIS function in any Python 3 library, so I do not need to upload to the DDBB all the data.
I tried with Shaperly but I can not find the way for the 3D distance


texts: How to use "geometryEngine" in arcGIS Version 3.9

texts:
We are using arcGIS javascript version 3.9 in our portal. Now we have new requirement to do reverse clipping on map. We have already written code to implement reverse clipping with the help of "esri/geometry/geometryEngine". The issue is that "esri/geometry/geometryEngine" was introduced in 3.13 and we are using 3.9 in our portal. When we try to load

require(["esri/geometry/geometryEngine"], function(geometryEngine) { /* code goes here */ });

in version 3.9, it gives error.
Is there any way to use "esri/geometry/geometryEngine" in version 3.9.


texts: Slow response of Open Data Cube for large amount of polygons

texts:
I have indexed 5 Sentinel-2 images and I am requesting the data (10 bands and the 5 images) for each one of 180,000 polygons included in a shapefile.
Each request needs 8 seconds for the related data to be retrieved. I use a 32GB VM with 8 cores.
Is there any alternative and more efficient way to make the request other than the code below ?
I tried also by loading firstly all the data and then slicing based on the polygon's bounds.
Finally, do you think multiprocessing will help?
shape_file = '/data2/cyprus/cyprus_parcels_2019_36SWD_3857.shp'
ds = fiona.open(shape_file)
crs = geometry.CRS(ds.crs_wkt)
for f in ds:
 feature_geom = f['geometry']
 geom = geometry.Geometry(geom=feature_geom,crs=crs)
 bounds = shape(feature_geom).bounds
 if 'MULTIPOLYGON' in geom.wkt:
    continue
 product = "s2a_sen2cor_granule"
 query = {
    'geopolygon':geom,
    'time': ('2019-01-01', '2019-12-31'),
    'product': product
 }
 data =
 dc.load(output_crs="EPSG:3857",measurements['B04_10m','B08_10m'],resolution=(-10,10),**query)


texts: Convert an object of class MCHu into a ShapeFile or to SpatialPolygonsDataFrame

texts:
I was creating Local Convex Hull Home Ranges in R with the adehabitat package and am unable to figure out how to convert them into Shape files to use in ArcMap.  The methods I've tried seem to require the object to be of class 'SpatialPolygonsDataFrame' or a couple others, but not class 'MCHu'(which is what the Convex Hull Home Ranges are).

How do I convert an 'MCHu' into a 'SpatialPolygonsDataFrame' or another method to convert 'MCHu' into a shapefile to export into ArcMap?


texts: Auto Incrementing Field based on Selection of Records in Feature Class?

texts:
I'm attempting to write a script that will add an incremented value to a field based on a selection of records.

For example, the script will select the first 300 records in a feature class and then using calculate field assign a value of 1 to a field called "Week_Assignment". It will then select the next 300 records and assign a value of 2 to the same field, etc.

I'm able to select the records, but not able to calculate the field with an incremented value. I found some generic code on an Esri help forum that will auto increment every row in a field, but I need to modify it to add the same value to all rows in a selection.

Here's the sample code:

rec=0
def autoIncrement():
 global rec
 pStart = 1
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec += pInterval
 return rec

autoIncrement()


texts:
I am trying to compare 2 fields(double) using python. I want to round to the nearest 10th and flag those with a "1" that are not equal. Here is what I have at the moment:

def Qcontrol(old, new):
  FieldA = math.round(old,1)
  FieldB = math.round(new,1)
  if FieldA == FieldB:
    return 0
  else:
    return 1

Return=
 Qcontrol(!old!, !new!)


texts:
I would like to host QGIS on a server so that people can browse my maps. Is this possible to do completely on a local network as some of the data is sensitive and cannot be uploaded online. Looking for something with the look and feel of ESRI's webapp builder via Portal/Enterprise.

Is this something that is possible/exists?


texts: Is Bundle adjustment a registration process?

texts: Since Image registration is the process of precisely overlaying two (or more) images of the same area through geometrically aligning common features (or control points) identified in the images, can the bundle adjustment be viewed as a registration process?


texts: QGIS Labeling automatic callout-line small polygons

texts:
I have a polygon layer and want to create a automatic labeling.
In this polygon layer are some pretty thin polygons. For this polygons the labeling should be placed outside of the polygon but with a callout-line points on the polygon.

How can I create such a labeling?


texts: How can I associate a consistent polygon area with a set of attributes using a topology?

texts:
Consider the following topology.

SELECT CreateTopology('topotestschema', 0, 0.01);
CREATE TABLE topotest (id SERIAL PRIMARY KEY, data TEXT);
SELECT AddTopoGeometryColumn('topotestschema', 'public', 'topotest', 'topogeom', 'polygon');
INSERT INTO topotest (data, topogeom) VALUES ('outer', toTopoGeom('POLYGON((0 0,0 10,10 10,10 0,0 0))'::geometry, 'topotestschema', 1, 0.01));

When we look at what's in the topology:

SELECT *, ST_AsText(topogeom::geometry) FROM topotest;

we get

╔════╦═══════╦═══════════╦═══════════════════════════════════════════╗
║ id ║ data  ║ topogeom  ║                 st_astext                 ║
╠════╬═══════╬═══════════╬═══════════════════════════════════════════╣
║  1 ║ outer ║ (3,1,1,3) ║ MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0))) ║
╚════╩═══════╩═══════════╩═══════════════════════════════════════════╝

But this gets changed around when I add more shapes:

INSERT INTO topotest (data, topogeom)
VALUES
('right', toTopoGeom('POLYGON((9 1,11 1,11 9,9 9,9 1))'::geometry, 'topotestschema', 1, 0.01))
,('inner', toTopoGeom('POLYGON((3 3,3 7,7 7,7 3,3 3))'::geometry, 'topotestschema', 1, 0.01))
;

Now the layer is:

╔════╦═══════╦═══════════╦══════════════════════════════════════════════════════╗
║ id ║ data  ║ topogeom  ║                      st_astext                       ║
╠════╬═══════╬═══════════╬══════════════════════════════════════════════════════╣
║  1 ║ outer ║ (3,1,1,3) ║ MULTIPOLYGON(((10 9,10 1,10 0,0 0,0 10,10 10,10 9))) ║
║  2 ║ right ║ (3,1,2,3) ║ MULTIPOLYGON(((10 9,11 9,11 1,10 1,9 1,9 9,10 9)))   ║
║  3 ║ inner ║ (3,1,3,3) ║ MULTIPOLYGON(((3 3,3 7,7 7,7 3,3 3)))                ║
╚════╩═══════╩═══════════╩══════════════════════════════════════════════════════╝

How can I keep it from modifying the polygon associated with existing rows when I add new ones?

Please feel free to say something if it seems like I'm on entirely the wrong track here. I'm basically experimenting with a topology and trying to figure out if it's useful for what I'm doing.


texts: What is the correct format for APA referencing maps produced in QGIS, including data layers?

texts:
I need to cite maps produced in QGIS. In addition to my own data, I have used the opensource MapQuest OSM basemap, and a 25m land type data layer obtained from Digimap.

I would also like to include the QGIS software itself. Is there a correct format for this? If not, what is the key information I need to include in the citation?


texts: Handling XML from within ArcGIS Server environment?

texts:
In my company, there is a web service that serves geographic data stored in a remote dababase through an URL. This service makes queries to the database and encodes the required data in XML format. So, to query data from a table (for example User), I do something like:

http://myCompanyName/service/QueryUserTable

So the service returns back to the client's browser a XML with all users stored in that table. Now, I would like to do the same from ArcGIS software: The thing would be connect this desktop application to that remote service, retrieve that XML file and parse it but all inside ArcGIS environment, not using other scripts or programs outside it. Is it possible? :S

Just in case it is not possible to do what I said, what is the best approach to handle and manage traffic in XML using ArcGIS?

It is a good idea and I like the way of retrieving the parameters like this long link you pasted. Didn't know about that module since I am pretty new with ESRI software.

But there is another constraint I did not mention in my first post. As the company is big, there is another "handicap" or problem: At this beginning stage I will have no control over the remote service, another developers group will "decide" how to encode the data. I mean, I think I cannot decide if the remote service is SOAP or not, REST or not and therefore, I don't know the XML tags they will use for the encoding. Maybe I can use the module you mention to connect to the service but I am not sure how data will be served on the other counterpart...

I guess I will be told the general schema of the files, and then I have to try to connect ArcGIS to that remote service and prepare some scripts to parse the incoming XML and allow non-IT users to use it in their layers...

Does it make sense? What do you think? It is a bit difficult to explain!


texts: Calculating maximum distance within polygon in x-direction (east-west direction) in PostGIS?

texts: I am interested in the maximum width of a polygon e.g. a lake, in east-west direction. Bounding boxes will help only in simple polygons but not in complex concave polygons.


texts: In QGIS, in the print composer screen if you save you work as 'Save as template' how do you reload it?

texts: I have various data that I am using in the print composer to create images and putting them into word. They will be marked soon by my assessor and I will likely go back and improve them, how do I re-access them?


texts: Managing the state of Python Add-In extension

texts:
I have a python add-in that includes a toolbar and an extension. A button on the toolbar should turn extension on and off.

A code for the button:

class Switch(object):
# Implementation for Sample_addin.btnSwitch (Button)
def __init__(self):
    self.enabled = True
    self.checked = False

def onClick(self):
    if not self.checked:
        self.checked = True
        print "Extension ON"
        # extExtension.__init__()
        extExtension.enabled = True
    else:
        self.checked = False
        print "Extension OFF"
        extExtension.enabled = False

A code for the extension:

class Extension(object):
# Implementation for Sample_addin.extExtension (Extension)
def __init__(self):
    self.enabled = True
    print 'Initialized'
# and so on...

And a part of config.xml:

<Extension  autoLoad="true"
            category="Tool"
            class="Extension"
            description="description"
            id="Sample_addin.extExtension"
            name="Extension"
            productName="Extension"
            showInExtensionDialog="true" />''

Everything seems to be OK. When I click the button, it freezes in 'checked' state and the extension gets marked in Extension dialog (Customize -> Extensions). But after this the extension doesn't start working (I have many functions in it like 'onCreate' and so on).

If I toggle it manually in Extensions dialog - everything is OK, it works.

And the reverse case - when the Extension works and I want to switch it off by clicking on the button, it only changes its state in Extension dialog but not in reality.

I've tried to reinitialize the extension by calling extExtension.__init__(). This gave me nothing (and no errors too, only this print statement).

How to track the state of an extension in python add-in and change it?


texts: How to display Google Maps as basemap in a Silverlight application?

texts:
I am writing a web mapping application and I use ArcGIS Server services. I want the application to use Google Maps as a base map. I am new to this job. I have searched a lot on the Internet. But I could not find anything useful. Is there anyone who could help me in this matter. Sample site are shown below.
http://esup.qom.ir


texts: Count points in polygon results in error "'NoneType' object has no attribute 'geometry' See log for more details"

texts:
I'm new to QGIS and have the following problem: I'm trying to use the "Count points in polygons" tool, but always get an error ('NoneType' object has no attribute 'geometry' See log for more details).

Uncaught error while executing algorithm
Traceback (most recent call last):
File
"C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py",line 203, in execute
self.processAlgorithm(progress)
File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\algs\qgis\PointsInPolygon.py", line 91, in processAlgorithm engine = QgsGeometry.createGeometryEngine(geom.geometry())
AttributeError:'NoneType' object has no attribute 'geometry'

If I use the same points layer, but a different polygon layer, everything works fine. So there seems to be a problem with my polygon layer.
Since I have no knowledge in python code, could someone please explain me, what I am doing wrong?


texts: arcpy.CalculateField_Management using sql statement

texts:
I am using python to calculate fields with a SQL statement. I need to select based on two attributes and calculate another attribute. The way I am typing it is incorrect syntax.

This is a small portion of the code and will run on a specific feature class.

Here is the code:

codeblock = """def calcFerry(FERRY_TYPE,AR_AUTO):
    if(FERRY_TYPE == 'B') and (AR_AUTO == 'N')
         return 'P'
    elif(FERRY_TYPE == 'B') and (AR_AUTO == 'Y')
         return 'A'"""

arcpy.CalculateField_management(outTemp_file, "FERRY",     calcFerry(!FERRY_TYPE!,!AR_AUTO!)", "PYTHON", codeblock)


texts: Split Polygon Layer by other Polygon layer in QGIS

texts:
I know there is a lot of previous posts to this subject, but I simply cannot find the solution I am after in the other posts.

I have to following problem: I have two very extensive polygone layers, one contains the planning zones of a region and one contains the actual allotments. Now one allotment can be in two differen planning zones, which is where I would like to split the allotment-polygons in to separet allotment-splitters per planning zone.

In theory, this sounds like a very basic thing, but I just dont seem to find anything in QGIS that enables it. I have thousands of polygons in both layers, so logically the Editing Tools where you have to select each pair of polygons/polygon-line is not an option. I wouldn't mind turning the planning zones into a multi-line-object and working from there but also that way I couldnt find a solution. Am I overlooking the obvious?

In ArcGIS this would be rather easy if one had the necessary licences, which I don't plus I am not a python-programmer.


texts: Geoserver and tomcat8 conflict of port numbers on cent os

texts:
I have installed geoserver binary and tomcat8 in /opt on centos 6.6.
My problem is, I can run geoserver only when I stop the tomcat service. How do I resolve the issue and start both on booting the system.


texts: Use new OpenLayers.Filter.Logical from OpenLayers into jQuery

texts:
I've a small GIS app, displaying stations on a map. Very simple.
At start, none of the stations are displayed, even if its layer is enabled.
What I want, when I click on a div, is to display one the station according to the OpenLayers.Filter.Logical  method. I plan to use jQuery to do that.
Problem : I don't master jQuery enough.
I'll provide you some codes I made. If you need mor to know, tell me.
My HTML code (extract)

 Projet

Station 1

JS code (extract)

station = new OpenLayers.Layer.Vector("station", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        protocol: new OpenLayers.Protocol.WFS({
            version: "1.0.0",
            url:  "http://585.182.555.82/sig/geoserver/wfs",
            featureType: "station",
            srsName:"EPSG:2154",
            featureNS: "http://www.obx.com/zodb",
            geometryName: "geom"
        }),

Example of a logical filter I already made directly into the station item (working) :

filter: new OpenLayers.Filter.Logical({
            type: OpenLayers.Filter.Logical.OR,
            filters: [
                new OpenLayers.Filter.Comparison({
                    type: OpenLayers.Filter.Comparison.EQUAL_TO,
                    property: "id_station",
                    value: "2"
                }),
]
        }),

You've understood it, my aim is to use the OpenLayers.Filter.Logical into the jQuery function.
I think I have to use something like this :
station.filter = new OpenLayers.Filter.Logical({
But I'm not sure how.

If I do

It's working (it's returning me "2", as according the logical filter)


texts: How do I query highways with overpass turbo?

texts:
I want to use overpass to show major highways for a US city. How can I query that from overpass turbo?

I tried the following code here and it returns tons of dots. I really want only lines showing the major highways

<!--
This has been generated by the overpass-turbo wizard.
The original search was:
“highway”
-->
<osm-script output="json" timeout="25">
  <!-- gather results -->
  <union>
    <!-- query part for: “highway” -->
    <query type="way">
      <has-kv k="highway"/>
      <bbox-query {{bbox}}/>
    </query>
  </union>
  <!-- print results -->
  <print mode="body"/>
  <recurse type="down"/>
  <print mode="skeleton" order="quadtile"/>
</osm-script>


texts: Converting DEM to point file gives ERROR 000870?

texts:
I have a DEM which has noisy parts to it. I am trying to (eventually) calculate the volume. However when I do cut/fill, there are obviously white sections where it hasn't been able to calculate the noisy parts. To get around this I was going to interpolate my DEM with one of the interpolate tools (IDW/Kriging/NN/Spline etc)

Problem: It has to be point file, so I was converting raster to point file when error 000870 appear stating that I already have an output file with the same name. I definitely don't and I have turned overwrite on just in case.

I was then going to convert back convert back to raster to do the volume.

Can anyone think of a simpler way or solve my mysterious file problem?

Using ArcGIS 10. All extensions possible, cannot download additional programs due to being on a university computer and not having admin rights.


texts: Does GPS measure faster speed better or lower speed better?

texts:
Assuming the same signal condition (HDOP, SNR, or other environmental factors), does GPS measure faster speed better? or lower speed better? Note that the speed measurement here is derived from Doppler effect, not from the position difference.
I think GPS measures high speed better (showing low relative error) because Doppler shift is too slight to be detected in low speed. I found a comment in the Stackoverflow post (texnic's comment in the AlexWien's post) questioning the same thing but it seems that there's no clear answer.
I tried to find research papers about relation between speed vs measurement accuracy of GPS but couldn't find anything handling this topic directly. So, can anyone give me the answer about it?
P.S.
This question was also uploaded in the stack overflow. Sorry if it violates the rule but I couldn't decided which one is the better place to post this question.
Edit: The "slower speed" means commonly < 10m/s


texts: Wrong object size after projecting spatial data from one coordinate system to another

texts:
I am having problems reprojecting items in ArcMap 10.1 from British National Grid To WGS84. (WKID: 27700 Authority: EPSG to WKID: 102100 Authority: ESRI)

All Polygon Layers I project [I use Project (Data Management)] from BNG to WGS are basically twice as big in metres as they were originally. Even the Shape Area and Shape Length inside the attribute table changes.

Any ideas on this would be greatly appreciated. Thank you.


texts: Resetting database coordinate reference system in QGIS Midvatten plugin?

texts:
When I started building a database in Midvatten I followed the tutorial:

https://github.com/jkall/qgis-midvatten-plugin/wiki/Tutorial

This tutorial suggests for setting up a new database:

"When prompted, select CRS 4326 (WGS84)"

As I was new to GIS I had not realised that this instruction was specific only to the demonstration data set that can be download to use with the tutorial and there might be reasons why I might choose another CRS.

I am now more up to speed with the question of Geographical versus Projected CRSs, (largely due to http://lyzidiamond.com/posts/4326-vs-3857) and now realise that I should have chosen EPSG 3857 rather than EPSG 4326.

The main reason is the warning in the Midvatten wiki for section plots:

"Please note! This feature assumes that projected coordinate reference systems are used (it will not work with tutorial data where CRS=EPSG 4326). It is also highly recommended to use the same Coordinate Reference System (CRS) for both the project and all involved layers."

I also now wish to use layers from the Openlayers plugin with my data set which also uses EPSG 3857.

I don't mind deleting and reloading the point observations data and stratigraphy data but as I have substantially changed the data domains it would take a significant amount of time to rebuild these.

Do I have to rebuild my database from scratch or is there a way of changing the CRS after the database has been created?


texts:
I am new to python ArcPy.
   I am trying to delete identical points in feature class.
the code is below:

>>> import arcpy
>>> fp=r"D:\NOCAS\TEMP\Belgum transistion.gdb\transistion\CenterlinePoints"
>>> sc=arcpy.da.SearchCursor(fp,["oid@","shape@"])
>>> dictionary={}
>>> for data in sc:
...     dictionary[data[0]]=data[1]
...
>>> uc=arcpy.da.UpdateCursor(fp,["oid@","shape@"])
>>> for row in uc:
...     for r in dictionary:

...         if row[0]!=r:

...             if sd[r].equals(row[1]):

...                 uc.deleteRow()

when trying to run the code it delete one point out five identical points and arcmap shutdowned.


texts: Installing ArcGIS Engine Runtime via image copy?

texts:
We have about 400 ArcGIS Runtime Engine 10.4.1 licences, where we will install on devices for emergency. These devices already up and running and using ArcGIS Engine Runtime 9.3. The time we are going to upgrade these device to the new version the time is critical for us to do that very fast. So we cannot just authorize each device separately.

Is there any way to do that via image copy?


texts: is there a way to reset the numbering in an attribute table in ARC 10.4? as in if features are deleted can the numbers be reordered in the new sequence left after deleting selections?


texts: Integrating ArcMap 10.1 Toolbox Model in R?

texts:
I have 30,000 raster files that needs to go through a small ArcMap 10.1 Toolbox Model I made (clip, flip, export to ASCII, etc.).  I have a R script which then imports all the exported raster file and does what I need it to do. Everything can simply be automated, except for the ArcMap 10.1 section which requires me to manually import everything via hand and run a Batch process (which frequently crashes the program and I lose the progress).

There has to be an easier way.  I looked into R-Toolbox-py, however it seems like it runs the entire R system inside ArcMap.  What I'm looking for is a simple way to process the required features (clip, flip, etc.) through the ArcMap Toolbox model. R-Bridge-Install tool by ArcGIS also doesn't work as I do not have access to 10.3.1.

Is there an easier solution to this issue without doing it by hand?  I can code in R, however I'm not very proficient in python.

ArcMap plays a very small role in my scripting/needs, therefore if there was an easy way to call the ArcMap Toolbox Model from R, that would work well because the toolbox model handles the input, and the output information without needing to even really open up the files.


texts: What is the difference between the the_geom and the_geom_webmercator in CARTO?

texts:
I'm trying to use CARTO for a school project and so I've just been trying to get the hang of the basics of using it and Postgis at the same time. I've read their docs about the difference between the_geom and the_geom_webmercator but I keep constantly tripping on it when writing queries.

Their docs:

  the_geom: field where the geometries (points, lines, or polygons)
  are stored in the CARTO dataset. The geometries must be in the
  projection EPSG:4326 (WGS84 - Geodetic coordinate system for World)

  the_geom_webmercator: this field is hidden in the dataset, and it
  stores the geometries in the projection EPSG:3857 (WGS84 Web Mercator
  (Auxiliary Sphere) - Projected coordinate system for World), this
  field is used to render the geometries on the map. This field is
  automatically updated when a geometry in the_geom field is added or
  updated.

Why does it seem like I find some query examples that use the_geom but others that use the_geom_webmercator? Why do I find queries that mix the two? Why do they have both of these to begin with? When using postgis functions like st_intersects, should I always call it on one or the other?

I've read up on projections and coordinate systems and I generally get all this. But for some reason, I am constantly confused when trying to write queries in the Carto context.


texts:
could someone please help with me with the following? Just need to get an idea of the flow of the code.

i) How do you use the ListFeatureClasses function to create a list of all the shapefiles. Loop through that list and use a copy tool to copy the datasets to another location. When this script is completed incorporates it as a script tool with parameters for the input and output workspaces?

ii) How do you modify the above script to incorporate an If statement. You only want to the first four shapefiles copied so once this amount of files has been copied use an if statement to halt to loop?

iii) How do you use the describe function and halt the execution of the program if an input layer is not of point geometry?


texts: Avoiding tiling/mosaicing in resulting image when Pan-sharpening Gram-Schmidt?

texts:
I mastered to georeference some declassified Hexagon satellite imagery from the 60s and would like to use it to pan sharpen a current Sentinel satellite image (already done as RGB). It is very important for me to retain as much color information from the original (the Sentinel) image as possible.

The Gram-Schmidt seems to be the best method for this purpose as it is the only one retaining enough color information from the RGB image.

When using ArcMap's "Image Analysis" tool to pan sharpen, the image looks perfect. But once exported to GeoTiff steps/tiling/mosaicing become aparent in the exportet image (see attachment).

I understand that this cannot be entirely avoided due to the method used that retains as much color information as possible but is there any way to make those steps more smooth/less obvious?


texts: Comparing rasters using Calculate Values in ArcGIS ModelBuilder?

texts:
I have two rasters with the same extents, each with all pixel values of the same value. Raster A has pixel values of 6 and Raster B has pixel values of 3. I know ModelBuilder has the Calculate Values tool that can process rasters with Python.

What code could I use to call these rasters and compare them to produce values?

What I am asking is if it is possible to essentially extract a number from a raster. So I know I can use Raster Calculator to generate a map where every pixel has a value equal to Raster A/Raster B. But I want to generate the number = to Raster A/Raster B so I can use it as an input in Modelbuilder.


texts: Snap Geometry to Layer (Line to Point) does not snap all lines within the tolerance limit

texts:
I've been using QGIS 3.14 version to snap geometries to layer.

When I try to run this, I have to fix the geometry first before I run select by location.
I checked the output, most are correct but there are lines that did not snap correctly within the tolerance limit I set like this: (The distance in the photo is not even 0.001)

Is the manual way the only way to fix this?


texts:
I did a topology check on a polygon layer. It turns out I have line errors "Must not have gaps". I don't understand as there are no gaps in the polygons.

I attached an image, the errors are the red lines. There is only one big polygon because I merged the islands with a Country.

Any idea on what could cause this error and how to fix it?


texts: Difference between layers using ArcGIS Geoprocessing?

texts:
I'm wondering if there is somewhere to get an output shapefile based on the difference between two layers. This would basically be the opposite of the clipping tool. for example. I have a roads polyline layer and a conservation easement polygon layer. I want to generate a new layer for roads that only includes those roads that are outside the easement polygon.

Is there some way to do this in ArcGIS Desktop?


texts: How to simulate way_pixels rules in GeoServer like in Mapnik?

texts:
Openstreetmap-carto project.mml have a queries like this:

way_area/NULLIF(POW(!scale_denominator!*0.001*0.28,2),0) AS way_pixels

Do I understand correctly that way_pixels just shows how many times the way_area is greater/less than the SD^2 area? And in code below from CartoCSS styles just mean "render marker if geometry area of a feature less than the SD area by 750 times"?

[feature = 'amenity_parking'] {
    [zoom >= 14][way_pixels > 750] {
        marker-file: url('symbols/amenity/parking.svg')
    }
}

Can I repeat the same behaviour in the GeoServer?


texts: Batch frequency of point groups to single table in ArcGIS Desktop?

texts:
I'm new in the GIS world.

I have a bunch of points (around 40,000), which are split in categories (1 to 9) and have a field which specifies their membership to a polygon (~50) (e.g, 30, which would mean that the point is located in polygon 30).

Some zones do not contain points of all categories.

I wish to know the frequency of each category per zone, in such a way that in the end I would have a table with the first column for the categories, and the zones listed in the first row.

Did somebody have the same challenge? How did you solve it?


texts: Can 3-opt local search algorithm be used for solving TSP based on network Dijkstra distances?

texts:
3-opt local search algorithm is widely used to solve Travelling Salesman Problem (TSP), as it provides near optimum solution in very less time as compared to Genetic Algorithm. What I understood that it works by swapping 3 edges at a time and swapping continues until the optimum solution is obtained. I don't understand how the initial swapping starts and what is the selection criteria for selecting the edges to swap for the next swap.
I read somewhere that it selects the edges that crossover each other and swap so that it does not.
Is it one of the method or only method for selection the edges for swap?
If that is not the only method, then what if I pass the actual road distances based on Dijkstra algorithm to distance matrix rather than Euclidean distances.
Would it be valid to use?


texts: Python GDAL: How to create virtual raster and cut out the polygon?

texts:
I would like to create virtual mosaic by gdal and then cut out the polygon and save it. The problem is I can't create virtual mosaic. How to solve that?
Error:
TypeError: invalid path or file: <osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x0000022773083B40> >

The code:
from osgeo import gdal
import geopandas as gpd
import shapely
import rasterio
import pandas as pd

df = pd.read_excel('C:/Users/ja/Inne/Desktop/warstwa_testowa.xlsx')

a = 'C:/Users/ja/Inne/Desktop/73776_931613_M-34-46-C-d-3-3.tif'
b = 'C:/Users/ja/Inne/Desktop/73776_933044_M-34-58-A-b-1-1.tif'
lista_tiff = [a, b]

completed_tif = gdal.BuildVRT('Vector_district.vrt', lista_tiff)

for index, row in df.iterrows():
    coo = shapely.wkt.loads(row['geometry'])
    multiplier = int(row['liczba_kon'])
    coo_buffer = coo.buffer(2 * multiplier)
    lokalnyid = row['lokalnyid']
    with rasterio.open(completed_tif) as src:
        out_image, out_transform = rasterio.mask.mask(src, [coo_buffer], crop=True)
        out_meta = src.meta
        out_meta.update(
            {
                "driver": "GTiff",
                "height": out_image.shape[1],
                "width": out_image.shape[2],
                "transform": out_transform,
            }
        )
        with rasterio.open(f'C:/Users/ja/Inne/Desktop/{lokalnyid}.tif', "w", **out_meta) as dest:
            dest.write(out_image)


texts:
I'm working on QGIS 2.14.
I'm using grass7:neighbors algorithm but i get artefacts on some sides of my images

I would like to clip my new image with the shape of the original to remove this problem, but as I use Vector -> Research -> Polygon from Layer I'm getting the boundaries as a rectangle like this

I succeeded getting the boundaries by reclassifying pixels values (outside raster values seems to have value 1) and then vectorize the raster. I tried to cut my first raster with those boundaries, but it's not working and I still get the rectangle when i try Vector -> Research -> Polygon from Layer.

When I try Raster->Convert and set the nodata value to 1 I get this

When I try the saga reclassification, and set values from 0.9 to 1 = 0/0 , I just get the rectangle filled with black.

Is there a way to tatally remove those values and to make them as the rest of the canvas ? (I don't want to give them another value because it would interfere with the grass:neighbors algorithm). Is ther a place where I can find what is the nodata value of my raster?


texts: Suppress legend item scaling (layers panel QGIS)

texts:
Is it possible to suppress the legend item scaling in the layers panel of QGIS 3.4.x?

I'ts a nice functionality, but confusing when using the magnifier > 100%


texts: Set the label size to fit the scale in the composer window

texts:
Is there a way to set the labels size of font to fit the scale of the map in the composer window but also to fit in the main window?

The map in the composer is size A0 and the scale is 1:5000.
In order to see the labels clearly in the composer, I set the labels font size to 23. So in the main window it's too large but in the composer it's OK.


texts: CreateFeatureclass management tool doesn't create Annotation Feature Classes?

texts:
I want to unify my gdb's, and later the layers in the mxd. I ran the following script successfully, but some of my layers are not the same as the template ones. Annotation Feature classes are not created. Why doesn't this tool create Annotation Feature Classes?

import arcpy, os, sys
#list of my desired FCs in a gdb
fcList = ['FC1', 'FC2', 'FC3']
gdb = "D:\\scripttest\\fillup\\test.gdb"
spatial_reference = arcpy.Describe("D:\\scripttest\\fillup\\test.gdb\\FC").spatialReference
has_m = "DISABLED"
has_z = "DISABLED"
workspace = "D:\\scripttest\\fillup"

for fcName in fcList:
    arcpy.env.workspace = gdb
    if not arcpy.Exists(fcName):
    #Searching for a template feature class in a folder to copy it's structure and appending it to fc[]
        fc = []
        while len(fc)<1:
            for dirpath, dirnames, filenames in arcpy.da.Walk(workspace, datatype="FeatureClass"):
                for fname in filenames:
                    if fname == fcName:
                        fc.append(os.path.join(dirpath, fcName))
        print fcName + " creating"
        template = fc[0]
        print ("temp: "+ template)
        #getting the geometry_type from the template
        desc = arcpy.Describe(template)
        geometry_type = desc.shapeType
        arcpy.CreateFeatureclass_management(gdb, fcName, geometry_type, template, has_m, has_z, spatial_reference)
        print (fcName + " created")
        del fc[:]


texts: Splitting feature using gaps around in ArcGIS for Desktop?

texts:
I have a shapefile which has only one feature. I want to separate it into a set of features. these features must be separated using gaps between each feature. How can I do it?

this is its image

I have used some tools such as erase and diffrence but they don't work because there is only one feature, So I must find a way to split this one feature using gaps into set of feature classes.

I am using ArcGIS 10.3 for Desktop.


texts:
For the record:  QGIS 3.2 and Mac OS 10.13.5

I am trying to load the symbols from https://www.mapbox.com/maki-icons/

I have tried creating /Users/rful011/Library/Application Support/QGIS/QGIS3/profiles/default/svg (Which is the default path for sgv files that I copied from Preferences->Option->System) and copying the .svg files into the directory but I still don't see any more symbols in the style dialogue.    (edit: because I was looking in the wrong place).

Where should the files appear in the style dialogue?

I then tried "open library" button -> import/export -> import ... and navigated my way to the directory where the svg files are:

As one can see from the screenshot the files are not selectable in the finder window.

After a bit more research (reading the qgis3 symbol library doc) I tried renaming some of the svg files to .xml and they then became selectable in the finder but when I tried importing them I got a dialogue box with the message "An error occurred during import: Incorrect root tag in style: svg".

I am starting to conclude that QGIS is not expecting one to import .svg files this way. ( edit:  this is definitely not the way to do it.  THis is for simple markers)


texts: CRS conversion using Lambert Azimuthal Equal Area project goes wrong for LAEA EUROPE using C#

texts:
When I try convert Lat Long from shapefile using Dotspatial for Europe LAEA and it gives null values for whole Europe region. It is a custom CRS provided by client

I have used Dotspatial and tried conversion based on given PROJCS and I have tried for Russia, North and South America  CRS and it works fine except Europe.

Below is the custom CRS provided by client

  PROJCS["WGS84LAEAEurope_54386_PCRS",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["FALSE_EASTING",4321000.0],PARAMETER["FALSE_NORTHING",3210000.0],PARAMETER["CENTRAL_MERIDIAN",9.999999999999989],PARAMETER["LATITUDE_OF_ORIGIN",52.0],UNIT["Meter",1.0],AUTHORITY["XXX",722486]]

It provided NaN for complete Europe region, like UK , Spain, Italy, France etc.,

Below are Lt and Long information which I try to convert

Lat= 57.04632319800004 Long =6.977629678000028
Lat= 57.04357274600005 Long =6.945070506000036
Lat= 57.04062228300006 Long =6.912269669000068


texts: I'm having trouble importing coordinates. I have a map over Scandinavia, and I need to import sampling coordinates that I've downloaded. The CRS I'm using is SWEREF99 TM EPSG:3006, and the format I'm importing is in SWEREF99 TM. For one of the files I get only one point, and for the other file I get absolutely nothing, no matter what I try. One of the other files I've got has the coordinates in lon lat, which I didn't have any trouble importing. Is there a way to easily convert the coordinates I have into decimal point of some sort? Or am I missing something?


texts:
I've got a polygon of agriculture field. I want to fill this polygon area with polylines in equal spaces between them (it supposed to represent the tree lines). those lines are not only for symbolic use, its should be for analysis.

Is there any way to convert this polygon into polyline like I described? I have tried to use the "polygons to line" tool but it gives me only the scope of the polygon and not the fill.


texts: How to add shapefile name as attribute during merge?

texts:
Possible Duplicate:
  How to add field with filename when merging shapefiles with ogr2ogr?

I have more 200 shapefiles with a unique name. I would like to merge all of these together and save the name of each individual shapefile as an attribute for each feature in the merged shapefile. Does anyone know if this is possible using the Open GIS (no proprietary GIS), plugins for Open GIS? Thanks!


texts: Visibility of SQL SERVER spatial data in QGIS depends on srid

texts: I loaded spatial data in ms sqlserver. Most data are in WGS 84 coordinates, but some are in Pulkovo 42. WGS 84 objects are visible well in QGIS. I checked it in 2.4 and 2.16 versions. But Pulkovo 42 features are invisible. It depends only on value of [geometry_columns].[srid] field. When I change value from 4284 to 4326, features became visible. Records in [spatial_ref_sys] don't affect the result.


texts: Setting Two Definition Queries as Parameters in Python Script tool causes: RuntimeError: LayerObject: Set attribute definitionQuery does not exist

texts:
I created a script that adds definition queries to every layer in an mxd depending on the layer name.

This worked fine as a stand alone script, but when I try to use it as a script tool with arcpy.GetParameterAsText() it fails on the second query.

Here's the snippet:

import arcpy

arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument(arcpy.GetParameterAsText(0))

def applyDefQuery(mxd):
    df = arcpy.mapping.ListDataFrames(mxd, "")[0]
    layers = arcpy.mapping.ListLayers(mxd, "", df)

    for layer in layers:
        if layer.name == "PROJECT_POLYGON":
            layer.definitionQuery = arcpy.GetParameterAsText(1)
        else:
            layer.definitionQuery = arcpy.GetParameterAsText(2)
    mxd.save()

applyDefQuery(mxd)

The idea is to prep an mxd for QA/QC when a subcontractor returns a checkout with a ton of layers. The script should select all of the project area polygons that they edited within and for every other layer select the records that the subcontractor edited.

I've tried troubleshooting this in a variety of ways: adding a second function for the second query, taking it out of the for loop and running it before the conditional statement, searching for null values in the editor fields instead of the ones I'm looking for (I was concerned that the value format - 'COMPANY\editorname' was causing an issue since it contained an escape character). Nothing seems to work except for leaving the third field in the GUI blank and just running the first definition query.

Here are the parameters for the script tool:

and the input I'm passing to it:

and the error:

Executing: QCPrep "C:\Users\bhamson\Documents\Company Check-    INs\VisualQC_04142015\VisualQC_04162015_2 (1).mxd" ""CONTRACT'" = 'AAA-CO-11-001_005'"     ""EDITOR" LIKE 'COMPANY%'"
Start Time: Tue May 05 18:30:38 2015
Running script QCPrep...
Failed script QCPrep...

Traceback (most recent call last):
  File "C:\Users\bhamson\Documents\Company Check-INs\VisualQC_04142015\QCPrep.py", line  17, in <module>
    applyDefQuery(mxd)
  File "C:\Users\bhamson\Documents\Company Check-INs\VisualQC_04142015\QCPrep.py", line  14, in applyDefQuery
    layer.definitionQuery = arcpy.GetParameterAsText(2)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\arcobjects\_base.py",  line 89, in _set
    return setattr(self._arc_object, attr_name, cval(val))
RuntimeError: LayerObject: Set attribute definitionQuery does not exist

Failed to execute (QCPrep).
Failed at Tue May 05 18:30:56 2015 (Elapsed Time: 17.58 seconds)

Any idea what the problem could be?


texts: How can I get a list of tables in my CartoDB account by size?

texts: The CartoDB dashboard only sorts tables by date modified or created. I would like a list of tables sorted by size.


texts: Shapefile creating "leaky" overlay in Google My Map after QGIS KML export

texts:
California has published, as shapefiles, maps of the various Fire Hazard Severity Zones in the state. I'm trying to use QGIS to convert these into a KML that I can then use as a Layer in a Google My Map. This is generally working well, however, on some Google pages and mobile apps, it seems like one of the KMLs is "leaking" or "bleeding" and showing a colored overlay on the entire map, and not just the bounded regions.

How can I clean up the file and properly view the Layer?

There are screenshots below. The steps to reproduce the issue are:

Download the State Responsibility Area (SRA) map shapefiles for Contra Costa county from the Cal Fire website.
Drag the entire fhszs7sn.zip file into the Layers panel a new project in QGIS (I'm using 3.8.2-Zanzibar), accepting the default CRS import settings.
Right click the new layer and go to Export -> Save Feature As. Select KML from the dropdown at the top, and leave the other settings as their defaults. (The KML file I'm having trouble with is here.)
Create a new Google My Map using the "+" icon at the lower right of this page. (My test map is here.)
Import the .kml file into a new Layer, click individual styles, and select HAZ_CLASS and select colors that are easy to distinguish.

At this point, everything should look normal. However, if you view the map in either/both of the following two ways:

Open the main Google Maps page, not the My Maps page above.
Click the three-line "hamburger icon" in the upper left, select "Your Places", then "Maps", then your just-created map. (I can't link to this directly for my sample map.)

Or:

Open the Google Maps app for iOS on your iPhone or iPad. (This is my primary use case.)
Tap the three-line "hamburger icon" in the upper left, select "Your Places", scroll the list at the top to the right if necessary, then tap "Maps", then your just-created map.

You should see the overlay covering the entire map. I'm able to distinguish by color that it's the "Very High" FSHZ that is causing the issue. I've tried using some QGIS tools to fix it, including "Check validity" and "Fix geometries", but neither seemed to help. What should I try?

I'm relatively new to GIS, so I apologize in advance if I've made a simply mistake. For what it's worth, the Alameda county LRA shapefile also had an error, but I was able to get QGIS to fix that one.


texts: How to remove a hole from a polygon in OpenLayers?

texts:
Is it possible to remove a hole from a polygon by using OpenLayers?
How can I do so?


texts: How to populate geography point using ST_MakePoint(long, lat)

texts:
I have a table called customer_campaign, and there is an empty column geom(geography) in it. So I want to fill it up by using ST_MakePoint(long, lat) in PostGIS, lat and long are two columns in the customer_campaign table. I tried the following query:

SELECT ST_MakePoint(customer_campaign.long, customer_campaign.lat);

But returns an error. How do I populate geography column with ST_Make(long, lat) in PostGIS ?


texts: QGIS - Automatically find and merge adjacent lines in large data set

texts:
I need to merge adjacent lines in QGIS 3.8, but am having trouble coming up with a way to identify and then merge them in a large data set. (3.8 is as far as my Mac can support)

I tried to convert my lines to polygons and use the "Join multiple lines" plugin. With both of these it seems like I'll need to select everything manually (not a feasible option) first or else it gets messy.

Ex. I want to locate lines like on the right here with touching endpoints, and merge them into one like on the left.

TIA!


texts: Difference between network dataset and geometric network in ArcGIS Desktop?

texts:
In ArcGIS Desktop, what is the difference between a "network dataset" and a "geometric network"?

Are there differences with what can be done with each?


texts: Drawing aerial routes for flights on Google Maps API between two markers?

texts: I have been drawing route till now using DirectionService. But now my current requirement is to draw a route on air.I did some research and found this link which helps us to draw curved polylines between two markers. But how do we make sure that this can be the exact flight route between two markers?. Completely Blank. What i want to do is something like this...Is there any good example to start with for drawing routes between two airports and animating a flight between the markers...Any help will be very much appreciated. Thanks in advance.


texts: PYQGIS: Summarizing values by defined field and write the result to a specific attribute

texts:
I'm using QGIS 3.10 with Python 3.7.
I would like to sum/group the length of 5 attributes by an ID and write the result in an other layer (to an attribute with the same ID).
Example:
I have two layers: one layer with houses and one with lines. Every house has 5 lines. Both layers have a field with an ID, so I can assign the lines to the houses.
I would like to sum/group the length of the 5 lines by the field "id" and write the result in a field of the house-layer.
I have a very large data set, so I need a solution with the shortest possible processing time.
my code (doesn't work / too long processing time (the computer crashes)):
for a in house_layer.getFeatures():
    houseID_a = a['id']
    dist_5 = 0
    for b in line_layer.getFeatures():
        lineID_b = b['id']
        if lineID_b == houseID_a:
            dist = b['length']
            dist_5 += dist
        a['sum_dist_5'] = dist_5
        house_layer.updateFeature(a)


texts: Select/export point features where the value change in two fields within the layer is less than a 5% change

texts: I have a point layer that has two fields representing the distance (in feet) of a point feature to coastline (old distance vs. new distance). I want to crosscheck both fields to see the accuracy of my new distance calculations compared to the previous. Is there a way to select/export all point features have distance value calculations with a less than 5% change in value from the old distance calculation field compared to the new distance calculation?


texts: Filtering out clusters of features based on collective area threshold using QGIS/python?

texts:
I have a map I generated that roughly estimates timber industry land ownership (designated by the green on the map) in the counties surrounding the Salish Sea in the US state of Washington:

I want to remove parcels below a certain area threshold (e.g "only show parcels larger than 20 acres"). But the problem with this is that many of the larger timber tracts are broken up into a bunch of smaller adjacent real property parcels. Thus, a timber company might have a contiguous 20 acre parcel that is broken up into four 5-acre parcels, and a simple "features where area > 20 acres" search would miss this parcel.

So what I am trying to figure out how to do is show any parcels that are either (a) larger than 20 acres or (b) part of a connected cluster/group of parcels that collectively cover more than 20 acres.

For example, zooming in on part of the above map:

I want to filter out any parcels that are smaller than the reference parcel, unless they are part of a larger cluster of adjacent/connected parcels (such as the one I labeled with the blue arrows in the image), and generate a new layer with these smaller, disconnected parcels removed.

Is there an easy way to do this through QGIS?

If not, does anyone have any ideas for Python libraries or command-line utilities I could use to achieve this?


texts: Vectorizing large-scale binary layer in Google Earth Engine?

texts:
I would like to vectorize binary raster layers in Google Earth Engine in order to be able to use the geometry of the resulting FeatureCollection in other tasks (e.g. clipping data or similar things).

One example: Let's create a vectorized map indicating high mountain regions based on SRTM data:

var srtm = ee.Image('CGIAR/SRTM90_V4');
var highMountains = srtm.gt(2000);
var mask = highMountains.reduceToVectors();

Unfortunately, this does not work - you'll receive the error message

Layer 2: Layer error: Image.reduceToVectors: Too many pixels in the region. Found 62208000000, but only 10000000 allowed.

Is there a way to achieve my goal? Or am I doing something wrong?


texts: Why is my attributes table greyed out for raster layers?

texts: I am trying to view attribute data for a SRTM GeoTIFF in QGIS, but for some reason I can't view the data.


texts: Positions from GPS lanes inside polygon using PostGIS?

texts:
I have a lot of lanes like the orange lanes in the below picture (latitude and longitude values).

Now I need all values from each lane inside the red polygon. What is the easiest way to realize this using PostGIS? I am not a GIS expert ;)

I would like read/load the files with lat/lon values, define the corners of the polygon (lat/lon??) and the output are all positions from each lane inside the polygon.


texts: Changing the transparency of a raster layer with ArcObjects

texts:
I'm trying to do some simple things with the .net framework and ArcObjects, but have little experience with this type of development.  I would like to iterate through the layers in a map document and change the transparency of the raster layers who's names begin with the letter "g".  Here is what I have so far as part of a windows form.

private void button1_Click(object sender, EventArgs e)
{
    IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
    IMap map = mxdoc.FocusMap;
    for (int i = 0; i < (map.LayerCount - 1); i++)
    {
        ILayer currentLayer = map.get_Layer(i);
        if (currentLayer.Name[0].ToString() == "g")
        {
            //current layer transparency = 50%
        }
    }
}

the line that is commented out is where I need help.  I can't seem to find the right code to do what I want within the ILayer class.  Is this not what I should be using?


texts: Create and populate field with name of feature (Select by Location) ArcPy

texts:
I am staring at 5600 individual points and 60 polygons. I need to create a new field in the point fc called 'Basin Name' and populate the field with the name of the Basin that it lies in. I am trying to develop a reusable script that can be used weekly to accomplish newly inputed points. I have two features, SSO which are points and Basins which are polygons. I have tried to follow this thread (Assigning field value based on location without join) but I keep receiving an error on the SelectLayerByAttributes.

The error I receive is "Runtime error Traceback (most recent call last): File "", line 19, in  File "j:\apps\arcgis\desktop10.3\arcpy\arcpy\management.py", line 7221, in SelectLayerByAttribute raise e ExecuteError: ERROR 000358: Invalid expression Failed to execute (SelectLayerByAttribute)."

import arcpy

arcpy.env.workspace = "Y:/SSO_BASINS/"

# Set overwrite option
arcpy.env.overwriteOutput = True

# KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)
# Create FeatureLayers
arcpy.MakeFeatureLayer_management("Y:/SSO_BASINS/Basins.shp", "lyr_Basins")
arcpy.MakeFeatureLayer_management("Y:/SSO_BASINS/SSO.shp", "lyr_SSO")

# Add a "BASIN" field
arcpy.AddField_management("lyr_SSO", "BASIN", "TEXT")

# Create a search cursor for the states
rows = arcpy.SearchCursor("lyr_Basins")
for row in rows:
arcpy.SelectLayerByAttribute_management("lyr_Basins", "NEW_SELECTION",
"\"OBJECTID\"  " + str(row.getValue("OBJECTID")))
arcpy.SelectLayerByLocation_management("lyr_SSO", "INTERSECT", "lyr_Basins",
"", "NEW_SELECTION")
arcpy.CalculateField_management("lyr_SSO", "BASIN",
"'{0}'".format(str(row.getValue("BASIN"))), "PYTHON_9.3", "")


texts: ArcSDE - spatial index grid size for empty table

texts:
I'm trying to write the SQL to create a feature class and spatial index in an Oracle SDE geodatabase using ST_Geometry. This is how I'm making the table:

CREATE TABLE STREETS (ID NUMBER(11) NOT NULL PRIMARY KEY, SEG_ID NUMBER(11) NOT NULL, SHAPE SDE.ST_GEOMETRY NULL)

Next I need to make the index which will look something like:

CREATE INDEX STREETS_IDX ON STREETS(SHAPE)
INDEXTYPE IS sde.st_spatial_index
PARAMETERS('st_grids=0 st_srid=4326')

The problem is I'm not sure what would be an appropriate grid size. I found this article which suggested using "one grid level with a cell size three times the average feature extent size", but the table is empty so I don't see how I would calculate that. Secondly, the spatial reference is WGS84 so what units would I be passing in? Degrees?


texts: Adding .tif layer to .mxd using ArcPy and showing RGB composite rather than stretched?

texts:
There seems to be no way to change the renderer using ArcPy.

Is there a workaround?

I am trying to show the RGB composite:

tif="test.tif"
newlayer=arcpy.mapping.Layer(rootDir+"\\"+tif)
# here is try to import a layer symbology to use, however it appears to have no effect on the output.
arcpy.ApplySymbologyFromLayer_management(newlayer, "C:\\mgmt\\tif_template.lyr")

arcpy.mapping.AddLayer(df, newlayer, "AUTO_ARRANGE")


texts: Custom Home Button In ArcGIS API for Javascript

texts:
Can you please take a look at This Demo and let me know how I can set up a custom home button

    var map;
    require(["esri/map", "esri/geometry/Point",
    "esri/symbols/SimpleMarkerSymbol", "esri/symbols/SimpleLineSymbol",
    "esri/graphic", "dojo/dom","dojo/on", "dojo/domReady!"], function (Map, Point,
    SimpleMarkerSymbol, SimpleLineSymbol,
    Graphic, dom, on) {
            map = new Map("map", {
            center: [0, 0],
            zoom: 5,
            basemap: "topo",
            slider: false, // set to false to remove default zoom buttons
            spatialReference: { wkid: 4326 }
        });
        on(dom.byId("zoomInBtn"), "click", function(evt){map.setZoom(map.getZoom()+1);});
        on(dom.byId("zoomOutBtn"), "click", function(evt){map.setZoom(map.getZoom()-1);});
        on(dom.byId("goHome"), "click", function(evt){});
    });


texts: Counting number of points in a point shapefile - ArcGIS 10.0 / arcpy

texts:
I assumed this would be an easy thing to do. I'm loading a point shapefile into a map like this:

layer = arcpy.mapping.Layer(file_name)

Is there a way to get the number of points in the shapefile? I see using the 'Describe' command that I can get some info about things... For instance,

arcpy.Describe(layer).table.fields

Gives me some info about the fields in the underlying table but I can't count the rows.


texts: Selecting polygons that are at least half covered by another layer using ArcGIS Desktop?

texts:
I have 2 polygon datasets. One is county-owned parcels within a county. The other is a land cover dataset that I have exported to only contain the areas that are forested. I would like to select from the parcel layer those parcels that are at least 50% covered by the forested areas polygon.

How could I achieve this using ArcGIS Desktop?


texts: Is there a way to create a reverse viewshed with ArcMap?

texts:
I am working on an environmental assessment. A data request has been made which it seems is kind of like a reverse viewshed, in that they want a layer that displays all locations within the raster that can view a single high point.

I have the appropriate raster and point layers, but if I run the Viewshed tool as is, pretty much the entire output is returned as not visible. Any thoughts on a workflow to try to capture what they are looking for? The only thing I've come up with so far is to use a bunch of random points surrounding the high point within the raster and include them in a viewshed, but this would not be a complete analysis of the area.

In this scenario, every cell in the raster would be a potential observer. I'll need the output to display values for those cells that can view the high point (1) versus those cells that cannot (0).

I tried running the standard viewshed tool, using the high point as the observer, along with the hillshade as the raster. The output returned the entire area as (0) with only a few cells directly around the high point as (1).


texts: RGEE getting pixel values and latitude and longitude

texts:
I am using Google Earth Engine on R. I'm trying to get a data set with the latitude and longitude of each pixel along with the value from the GEE Image Collection.
I know its been done in JavaScript - I just can't translate it over.
library(scales)
library(sp)
library(maptools)
library(rgdal)
library(cshapes)
library(rgee)
  ee_Initialize()
  #ee_reattach() # reattach ee as a reserve word
  #Adds a band containing image date as years since 1991.

 # Map the time band creation helper over the collection.

  collection <- ee$
    ImageCollection('JRC/GHSL/P2016/SMOD_POP_GLOBE_V1')$
    select('smod_code')

  collection2 <- collection$filter(ee$Filter$date('2010-01-01', '2015-02-01'))

##### Creating the Shapefile
b <- cshp(date = as.Date("1990-01-01"))

SouthAfrica <- c()
SouthAfrica2 <- c()
SouthAfrica<- b@polygons[[154]]@Polygons[[2]]@coords
  for(j in 1:nrow(SouthAfrica)){
    SouthAfrica2 <- c(SouthAfrica2, SouthAfrica[j,])
  }
SouthAfricaShap <- ee$Geometry$Polygon(SouthAfrica2)

##### Beginning the Conversion #####

SouthAfrica2000 <- collection2$filterBounds(SouthAfricaShap)
ImageSouthAfrica2000 <- SouthAfrica2000$first()
ImageSouthAfrica2000 <- ImageSouthAfrica2000$select('smod_code')

coordsImage = ee$Image$pixelLonLat()$reproject(ImageSouthAfrica2000$projection())
joinedImage = coordsImage$addBands(ImageSouthAfrica2000)

valuesList <- joinedImage$reduceRegion(c(
  reducer= ee$Reducer$toList() ,
  geometry=  SouthAfricaShap
))$values()$get(0)

print(valuesList)

The value list function is what's not working! I would love any help I can get!
THE ERROR THAT SHOWS UP IS REALLY MESSY AND SAYS THIS:
Error in py_call_impl(callable, dots$args, dots$keywords) :
  EEException: Invalid argument for ee.Reducer(): ({'reducer': <ee.Reducer object at 0x1241129e8>, 'geometry': ee.Geometry({
  "functionInvocationValue": {
    "functionName": "GeometryConstructors.Polygon",
    "arguments": {
      "coordinates": {
        "constantValue": [
          [
            [
              13.20889117993994,
              -16.97360893198885
            ],
            [
              13.260827179988269,
              -16.984863931999314
            ],
            [
              13.317782180041377,
              -16.97805493199303
            ],
            [
              13.34630918006792,
              -16.970672931986087
            ],
            [
              13.376036180095525,
              -16.970908931986344
            ],
            [
              13.472082180185055,
              -17.010835932023554
            ],
            [
              13.494300180205641,
              -17.02555493203721
            ],
            [
              1


texts: Getting tiles from WMTS document using GeoTools?

texts: I was wondering if there is a library which I can use to retrieve TileMatrix/Tile information from a WMTS Capabilities document. I have been trying to use geotools but it seems as though it does not support WMTS (just WMS).


texts:
using QGIS 2.2 and trying to figure out a way to obtain a distance measure for the width of pavement.
I have 2 shapefiles: one of points which are at 5m intervals along the pavement centreline and the second is the polyline of the pavement.
I'd like to take a distance measure from the point to both edges of the pavement, thus giving me the width. The real challenge is to get a width measurement that is the least angular change and the maximum distance - this is in order to deal with curved and wide pavements as seen in the graphic.
Advice will be most welcome - unfortunately, I can't code which is quite a limitation when trying to get to this solution.

Many thanks,
Khan


texts: Import a print template from AutoCad into QGIS Print Composer

texts:
I use QGIS (version 2.18.19) and my consulting colleagues use AutoCAD. A lot of our clients ask me to put my QGIS maps into AutoCAD format. Is there a way to simply import the print composer template from a CAD file directly into QGIS to then use?

Currently I have been saving as a PDF and then the CAD guys here import into CAD and create it that way. It is very time consuming and relies on other people.

If I could get the "look" of the AutoCAD format for QGIS then I could continue to happily make my beautiful maps in QGIS and not have to try and convert into AutoCAD.


texts: PyQGIS spatial operators(contains,equals,etc.) w.r.t. layer geometries(polygon,point,line,etc.)

texts:
I am trying to implement QGIS spatial query functionality in my plugin. I want to eliminate spatial operators which are not suitable for certain geometries. For eg. when I select layer with point geometry only certain(crosses,intersects, is disjoint,touches, Within) spatial operators are given by QGIS. Want to determine how to get to this result.

Got the layer geometry types using:

    layer = iface1.activeLayer()

    if layer.wkbType()==QGis.WKBPoint:
        print 'Layer is a point layer'

    if layer.wkbType()==QGis.WKBLineString:
        print 'Layer is a line layer'

    if layer.wkbType()==QGis.WKBPolygon:
        print 'Layer is a polygon layer'

    if layer.wkbType()==QGis.WKBMultiPolygon:
        print 'Layer is a multi-polygon layer'

    if layer.wkbType()==100:
        print 'Layer is a data-only layer'


texts: Qgis 1.8.0 throws error "Can not find the entry point sqlite3_rtree_geometry_callback"

texts:
Loading QGIS on Window 8 causes this problem:

Can not find the entry point sqlite3_rtree_geometry_callback in the procedure of dynamic loading C:\....\Quantum GIS Lisboa \ bin \ spatialite.dll

Everything worked properly for many months.


texts: Is it possible to disable an extension addin on start-up?

texts:
I am currently creating an addin, that can be toggled on/off from the extension manager. What I want to do is when ArcMap is being loaded, to check the license level being used (which I have done), what I then would like to happen is that the extension in the Manager is disabled and the text changed to something like "xxx ArcEditor Required".

Is this possible to do? I know I can change the name of the Extension via the config file, but can't seem to find the property to do it outside of it.


texts: Double or triple 45 degree hatch in GeoServer - possible?

texts:
My end goal is to have a style like this in GeoServer:

I have tried using ExternalGraphic in SLD but that only seems to work with vertical/horizontal hatch.

Edit:
By using ExternalGraphic with 45 degree hatch i get a offset like:

sld:

<sld:ExternalGraphic>
  <sld:OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="file:/C:/Users/asawea/Desktop/1012.png"/>
  <sld:Format>image/png</sld:Format>
</sld:ExternalGraphic>

Tips?


texts: QGIS Issue with 3D view altitude binding

texts:
I'm running into an issue with the 3D view in QGIS.

I have a DEM, which is displaying correctly and a building footprints polygon layer from OSM. It is draping correctly as a 2d polygon on the DEM in the 3D view, but when I go into the 3D View settings and enable 3D view rendering for the buildings layer to 'extrude' them, the 3D polygons are displaying far above the DEM vertically. I have switched the Altitude Binding to Terrain, which I believe should adhere them to the DEM since I have the DEM defined as the terrain in the 3D view, but nothing changes.


texts: Nesting project file name into layout in QGIS

texts:
Is there any way to automatically take information from the project file name into several boxes on a layout in QGIS?

For example, if I had a file name of "185.005R1 - Geology Map.qgz" how could I split it into separate textboxes, e.g. the Drawing number ("185.005"), the title ("Geology Map"), and the Revision number ("1")?


texts: OGR2OGR and SQL: selecting 30 biggest places from a country?

texts:
I use ogr2org to filter naturalearth places.shp :

ogr2ogr -f GeoJSON -where "ADM0NAME = 'India' AND POP_MAX > 500000" places.tmp.geo.json ne_10m_populated_places.shp

I would prefer to select the 30 biggest places (via ORDER BY POP_MAX DESC LIMIT 30). In the ORG SQL doc recommands something such :

SELECT * FROM property WHERE class_code = 7 ORDER BY prop_value DESC

I tried both

ogr2ogr -f GeoJSON -where "ADM0NAME = 'India' ORDER BY POP_MAX LIMIT 1,30" places.tmp.geo.json ne_10m_populated_places.shp
ogr2ogr -f GeoJSON -sql "SELECT * FROM ne_10m_populated_places WHERE ADM0NAME = 'India' ORDER BY POP_MAX LIMIT 30" places.tmp.geo.json ne_10m_populated_places.shp

But it doesn't work. What does I do wrong ? How to query right ?


texts:
Several county residents have been contacting our GIS office to find out why they are not receiving packages or why the delivery person gets lost finding their address.  I know it sounds like an odd question to as GIS, but it turns out thier packages are being delivered by UPS, FedEx, etc. All of them use Google Maps to locate an address in our county.  I checked a few of the addresses in question and they are either non-existent or in the wrong location.  I also found several subdivisions where not a single road is showing up in Google Maps either.  Now I know I can submit a "suggestion" to Google for them to update the map, but I really do not want to do this several hundred times to correct all the errors I have found so far.
Is there a way to upload GIS data from the county in order to update their maps? If so, is there an app, interface, or procedure available to do so?


texts: calculation precision during projection transformation

texts:
I am using wgs 84 base tiles, and the tile is not rendered as image but with geojson format.

I try to draw them in the tile according to the data. And the geometry of the returned feature is lng and lat format.

This is the formula I used:

double tileNumWidth = 1 << (zoom+1);
double tileNumHeight = 1 << zoom ;
float px = (float) ((tileNumWidth * fromLngToX(x) - tileX) * 256);
float py = (float) ((tileNumHeight * fromLatToY(y) - tileY) * 256);

public static double fromLngToX(double lng) {
    return (lng + 180) / 360d;
}

public static double fromLatToY(double lat) {
    return (90 - lat) / 180d;
}

It worked as expected in a lower zoom level:

http://pbrd.co/1qBajKH

However once the map is zoomed to a higher level like 18+, then I will found that the lines cross two tile may be not conjugated:

http://pbrd.co/1qBaN3k

You may notice the obvious drift(not much, seems only 1 px).

So I wonder if this is caused by the precision of calculation?


texts: Sublayers/subgroups under single item in ArcGIS Server

texts:
We would like to develop hierarchies for viewing layers in our ArcGIS Server based IMS. The best solution to do this has been so far to create groups and sub-groups in ArcGIS Desktop before publishing services. However, with this approach we cannot create groups under single items (shapefile/feature class/raster). We would like to do this as it would reduce the amount useless groups and thus useless headers. Is this possible and how could it be done?

I might also be on completely wrong tracks, when I try to do this in Desktop before publishing the service. Feel free to share if you have any other approaches.


texts: Calculate distance in km from point to polygon (given in lat/long) using ArcGIS?

texts:
I have two datasets, one of polygons and the other of points, both spread across the United States. Both have decimal degree coordinates in WGS 1984. I want to determine the distance from each polygon to the three nearest points in km. My initial approach was to use the Generate Nearest Table tool but because the inputs are in lat/long degrees, that is how the outputs are given.

This is a very similar question to one asked previously:
Calculate distance in km to nearest points (given in lat/long) using ArcGIS Desktop and/or R?
The difference is point to polygon vs. point to point. So I don't know how the code given in the answers would work. I'm not too familiar with R either so any desktop based solution would be more straightforward for me.

The Generate Nearest Table seems perfect for this it just doesn't work because of the lat/long coordinates I'm starting with. Should I reproject the datasets into something that would allow distance calculation in km?


texts: Setting feature attribute by name via QGIS python api?

texts:
I try to set feature attributes by attribute name with this code

pr = layer.dataProvider()
pr.addAttributes([QgsField("test", QVariant.Int)])
layer.updateFields()
for feature in layer.getFeatures():
 attrName = 'test'
 feature[attrName] = 1

but in result I have NULL in all fields.
If I use

feature.setAttributes([1])

it works properly.

What I do wrong? Why feature[attrName] = 1 doesn't work?

Update:
Find this solution

pr = layer.dataProvider()
pr.addAttributes([QgsField("test", QVariant.Int)])
layer.updateFields()
for feature in layer.getFeatures():
 attrName = 'test'
 feature[attrName] = 1
 pr.changeAttributeValues({feature.id() : {pr.fieldNameMap()[attrName] : 1}})


texts: How can I ensure the image exported from Print Composer is accurate?

texts: I'm using QGIS version 2.12.2, and using print composer to export my map to an image file (I've tried pdf as well). In the print composer it displays fine - just a basemap overlayed with land type data. However, as soon as I export the map, the basemap and the overlayed data are misaligned; rivers are categorised as fields etc. How can I ensure that the image produced from print composer is accurate?


texts: Google Earth Engine: batch image export with JavaScript API

texts:
My goal is to export a large mosaic into tiles defined by individual features in a vector grid (FeatureCollection).
GEE debugging guide states that

Export is a client-side function and won't work in mapped functions.

But when I use client-side loop the whole thing freezes for a long time due a to multiple .getInfo() requests:
// vector_grid is a FeatureCollection with 'tile_id' attribute

var tile_list = vector_grid.getInfo()["tile_id"];

for (var i in tile_list) {
  var tile = ee.Feature(tile_list[i]);
  var tile_id = feat.get("tile_id").getInfo().toString();
  var tile_geom = tile.geometry();

  Export.image.toDrive({
    image: mosaic,
    folder: 'export_path',
    fileNamePrefix: tile_id,
    region: tile_geom,
  });
}

The question is:
Is it possible to do a batch export defined by FeatureCollection in a more rational way within the JavaScript API?


texts: Concatenating text and string in Field Calculator expression (python)

texts:
I'm trying to copy an attribute value from one field to another (sum of raster values), but I'm having issues in writing a functioning expression for a Calculate Field script.  I'm adding a field called SUM and trying to copy the value from the field "SUM_" + filename[0:-4]. I have multiple files I am trying to do this on and the field names are different in each, which is why it is being written in this manner.

Here's the code I have so far:

for table in arcpy.ListTables():
    arcpy.AddField_management(table, "SUM", "DOUBLE")
    expression = ["SUM_" + str(filename[0:-4])]
    arcpy.CalculateField_management(table, "SUM", expression, "VB")

This doesn't work as I receive the following error:

File "C:\Users\kellyj\Desktop\Projects\CostAnalysis\Route_Cost_Analysis.py", line 82, in <module>
    arcpy.CalculateField_management(table, "SUM", expressionvalue, "VB")
File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 3354, in CalculateField
    raise e
RuntimeError: Object: Error in executing tool

Any thoughts as to how I can fix this?