texts:
Note: This question is specifically about installed, desktop software. There is another question specifically about free cloud-based software and services.

What free programs should every GIS user have installed?

I'm not necessarily referring to ESRI extensions or open-source products, but others that increase your productivity and ability to handle GIS tasks.

For example:

Notepad++ for writing code snippets or editing XMLs. Paint.NET or GIMP for quick graphic editing.
I use Google Tasks daily and I think it's worth mentioning. It's not GIS-specific, but it's a great tool, especially if used independently and on multiple projects where purchasing time-management software isn't reasonable.
While it's not focused on GIS development, Rainmeter has proven to be very useful in terms of increasing productivity and monitoring system resources. I have created a GIS "sidebar" on my desktop that holds all of my development tools, as well as links to the online resources I used the most. It's nice to be able to use one location, rather than many (e.g. taskbar, bookmarks in browser, search engine).


texts: Resampling using gdalwarp gives me no data

texts:
I have an image with the following information.

Driver: GTiff/GeoTIFF
    Files:test_DEM.tiff
    Size is 91, 73
    Coordinate System is:
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]]
    Origin = (-87.238000000000000,40.620599999999996)
    Pixel Size = (0.000100000000000,-0.000100000000000)
    Metadata:
      AREA_OR_POINT=Area
    Image Structure Metadata:
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  ( -87.2380000,  40.6206000) ( 87d14'16.80"W, 40d37'14.16"N)
    Lower Left  ( -87.2380000,  40.6133000) ( 87d14'16.80"W, 40d36'47.88"N)
    Upper Right ( -87.2289000,  40.6206000) ( 87d13'44.04"W, 40d37'14.16"N)
    Lower Right ( -87.2289000,  40.6133000) ( 87d13'44.04"W, 40d36'47.88"N)
    Center      ( -87.2334500,  40.6169500) ( 87d14' 0.42"W, 40d37' 1.02"N)
    Band 1 Block=91x22 Type=Float32, ColorInterp=Gray
      NoData Value=-3.40282001837565598e+38`enter code here`

I am trying to resample this image to a 5m x 5m resolution. I use the following command:

 gdalwarp -t_srs EPSG:4326 -tr 5 5 -tap -r bilinear -dstnodata -3.40282e+38 test_DEM.tiff DEM_result.tiff

But no image is being created. I am not sure where I am going wrong. This is the message I am getting:

Creating output file that is 1P x 1L.
Processing input file
test_DEM.tiff.
Using internal nodata values (eg. -3.40282e+38) for image test_DEM.tiff.
0...10...20...30...40...50...60...70...80...90...100 - done.


texts:
Does anyone know where I can get free land use data for Manhattan? I'm specifically looking for fraction of residential, commercial and other landuse data for each census tract for the New York County.

I did search the NYC - Bytes of Big Apple website and found that they charge an appalling 300$ per borough for the land use at tax lot level.

Any help is greatly appreciated...


texts: How to calculate Distance, Azimuth and Dip from two XYZ coordinates

texts:
I have two XYZ coordinates and I would like to calculate the Distance, Azimuth and Dip between them.

I have come across a similar question on this site that calculates it the other way around (with a given starting coordinate and a given Distance, Azimuth and Dip and it calculates an end coordinate.)

I have been trying to figure it out and I am sure its very simple.

Here is that other question page.
It says that this is a conversion between a spherical and Cartesian coordinates.

How to convert Distance, Azimuth, Dip to XYZ?


texts: How to determine nearest parcel neighbors, and adding direction info (N,S,E,W) in ArcPy?

texts:
I am working with a cadastral project, and I need a solution for a certain process, neighboring. Based on ArcGIS.

Scenario: Have some parcels, every parcel has some neighbors.

I would like to add info in every parcel for:

1-Nearest neighbor grouped by direction (N,S,E,W)

import arcpy

try:
    #Read parameters from dialog
    fc = arcpy.GetParameterAsText(0)
    fieldName = arcpy.GetParameterAsText(1)

    #Add appropriate fields if they don't already exist
    if fieldName+"_N" not in arcpy.ListFields(fc):
        arcpy.AddField_management(fc, fieldName+"_N", "TEXT")
    if fieldName+"_S" not in arcpy.ListFields(fc):
        arcpy.AddField_management(fc, fieldName+"_S", "TEXT")
    if fieldName+"_E" not in arcpy.ListFields(fc):
        arcpy.AddField_management(fc, fieldName+"_E", "TEXT")
    if fieldName+"_W" not in arcpy.ListFields(fc):
        arcpy.AddField_management(fc, fieldName+"_W", "TEXT")

    #Make a feature layer for SelectByLocation
    lyrFile = "lyr"
    arcpy.MakeFeatureLayer_management(fc, lyrFile)

    #Function that creates a point feature and performs neighbor selections and returns adjacent page values
    def findLocation(x, y):
        arcpy.CreateFeatureclass_management(r"in_memory", "tempFC", "POINT")
        ptLayerCur = arcpy.InsertCursor(r"in_memory\tempFC")
        ptLayerRow = ptLayerCur.newRow()
        pt = arcpy.CreateObject("Point")
        pt.X = x
        pt.Y = y
        ptLayerRow.Shape = pt
        ptLayerCur.insertRow(ptLayerRow)

        arcpy.SelectLayerByLocation_management(lyrFile, "INTERSECT", r"in_memory\tempFC", "", "NEW_SELECTION")

        lyrCursor = arcpy.SearchCursor(lyrFile)
        lyrRow = lyrCursor.next()
        pageNumber = ""
        while lyrRow:
            #pageNumber = lyrRow.MAPSHEETNUMBER
            pageNumber = getattr(lyrRow, fieldName)
            lyrRow = lyrCursor.next()

        arcpy.Delete_management(r"in_memory\tempFC")

        return pageNumber

    #Loop through each index polygon and calculate appropriate neighbors
    mapIndexCursor = arcpy.UpdateCursor(fc)
    mapIndexRow = mapIndexCursor.next()
    while mapIndexRow:
        #Calculate North Page XY values
        northX = (mapIndexRow.Shape.extent.XMin) + ((mapIndexRow.Shape.extent.XMax - mapIndexRow.Shape.extent.XMin) / 2)
        northY = (mapIndexRow.Shape.extent.YMax) + ((mapIndexRow.Shape.extent.YMax - mapIndexRow.Shape.extent.YMin) * 0.15)
        setattr(mapIndexRow, fieldName+"_N", findLocation(northX, northY))

        #Calculate South Page XY values
        southX = (mapIndexRow.Shape.extent.XMin) + ((mapIndexRow.Shape.extent.XMax - mapIndexRow.Shape.extent.XMin) / 2)
        southY = (mapIndexRow.Shape.extent.YMin) - ((mapIndexRow.Shape.extent.YMax - mapIndexRow.Shape.extent.YMin) * 0.15)
        setattr(mapIndexRow, fieldName+"_S", findLocation(southX, southY))

        #Calculate East Page XY values
        eastX = (mapIndexRow.Shape.extent.XMax) + ((mapIndexRow.Shape.extent.XMax - mapIndexRow.Shape.extent.XMin) * 0.15)
        eastY = (mapIndexRow.Shape.extent.YMin) + ((mapIndexRow.Shape.extent.YMax - mapIndexRow.Shape.extent.YMin) / 2)
        setattr(mapIndexRow, fieldName+"_E", findLocation(eastX, eastY))

        #Calculate West Page XY values
        westX = (mapIndexRow.Shape.extent.XMin) - ((mapIndexRow.Shape.extent.XMax - mapIndexRow.Shape.extent.XMin) * 0.15)
        westY = (mapIndexRow.Shape.extent.YMin) + ((mapIndexRow.Shape.extent.YMax - mapIndexRow.Shape.extent.YMin) / 2)
        setattr(mapIndexRow, fieldName+"_W", findLocation(westX, westY))

        mapIndexCursor.updateRow(mapIndexRow)
        mapIndexRow = mapIndexCursor.next()

except Exception, e:
  import traceback
  map(arcpy.AddError, traceback.format_exc().split("\n"))
  arcpy.AddError(str(e))
  arcpy.Delete_management(r"in_memory\tempFC")


texts: QUAC atmospheric correction for Landsat-8 (ENVI 5.0)

texts:
I try to realize QUAC atmospheric correction for the band 4 and band 5 of Landsat-8 in ENVI5.0. This version doesn’t support Landsat8 (doesn’t have parameters of Landsat-8).
At the first step I glued Band3, Band4 и Band5 by Layer stacking Tool in a single file.
Then (as I understood) I need to put some parameters in this file by dint of Edit Envi Header Tool.
I put Wavelenght from the official page:
http://landsat.usgs.gov/band_designations_landsat_satellites.php

The program doesn’t know the sensor type, so it remained “Unknown”.

I don’t know any other parameters for Landsat-8 and which parameters  I need to include.
So, I realized QUAC correction. But the result of NDVI calculation was failed.

Which parameters of Landsat-8 for QUAC atmospheric correction I need to include?


texts: Simple workflow for a GPX track export from GIS data, avoiding illogical geometry (on QGIS)?

texts:
GPX exports from GIS information are difficult, because a normal dataset - made of a number of separate line features - does not specify the order in which lines should be interpreted.

A good GPX track has an order to it. It's designed to represent an actual journey (direction, order, etc are inherently part of the format). It may have only one 'track' or separate track 'segments'. If it does have segments the expectation would (often) be that one segment follows from another in some logical spatial sense - perhaps the end of one segment is followed by the beginning of the next nearby segment implying a break in a journey, but not a teleportation to another location.

Exporting from GIS into a GPX track, at least when using QGIS, typically introduces unwanted track segments - and may not present these segments in the GPX file in a logical order. Devices interpreting a GPX file may be unforgiving, and avoiding such issues is important.

I often hit issues with GPX export. For example, to avoid creating separate track segments it's important to snap lines, then combine them - but issues like duplicate nodes cause problems. Line direction also matters. As does any crossing of lines (entirely acceptable in GPX - less so in GIS).

Is there a reliable workflow (perhaps involving other open source software too) which will reliably produce a tidy GPX track without illogical track segments? This should be able to cope with a starting point with standard GIS line data - multiple features- perhaps not correctly snapped at line ends.


texts: Can a routing service take polygon geometry (for stops)?

texts: I have built a small road featureclass, created a network off of it using network analyst, and created routes successfully in ArcMap.  I'm now trying to create these routes in a silverlight web application and in order to do that I need to pass the geometry of two stops (polygons) to the routing service.  I'm not having any luck and I think it's because the centroid being calculated by the web application isn't "exactly" on the routing service.  So I'm wondering instead of passing a point centroid, I could instead pass the polygons geometry of my two stops to the routing service?  Does anyone have an example of how to do this?  Or is this possible?


texts: PostGIS with GeoJSON: `geometry` v. `geography`

texts:
Does the GeoJSON spec imply or strongly suggest that PostGIS geometry type should be used over geography?

Looking at section 3.1.1 I am seeing that

  In other words, every point on a line that does not cross the
     antimeridian between a point (lon0, lat0) and (lon1, lat1) can be
     calculated as

  F(lon, lat) = (lon0 + (lon1 - lon0) * t, lat0 + (lat1 - lat0) * t)

  with t being a real number greater than or equal to 0 and smaller
     than or equal to 1.

This would seem to assume a flat/cartesian coordinate reference system? Otherwise, the points should lie along a great circle arc?


texts: OGR/Python- Check whether or not points are in a polygon

texts: I have two shapefiles, one points file and one polygon. They are in the same projection (Lat/Long). I would like to run a quick command-line command or a Python script with the OGR module that creates a new points shapefile of only the points that are within that polygon. Is there a good way to do this?


texts: Chronological sort & indexing problems with QGIS Time Manager

texts:
I am trying to use QGIS Time Manager to create an animated map of property connected to polygon-geometry. The one article on animating polygons I've found on this site hasn't helped because my data is joined to an existing shapefile, and the spatial index (FID) of the target file is preserved, meaning it won't resort chronologically to play nicely with the Time Manager plugin.

As a result, the Time Manager does something weird where it applies an auto-filter Query to my layer's data and arbitrary zero end date of '0001-01-01'. I did not set an "end date", but used "Accumulate features" setting

I have tried re-indexing using the Add Autoincremental Field to create a new index (with the sort condition set to Ascending on the Date field), but the FID of the new shapefile persists... Am I using the right strategy?


texts: Using ArcGIS Desktop with SQL Server Spatial as backend?

texts:
I'm wondering how best to use ArcGIS Desktop and SQL Server 2008 together. At the moment we run manual imports and exports between sql 2008 and arcgis 10 shape files. With sql 2008 having a geography data type I was expecting it to be simple to replace to the map files with a relational database, but for some reason it doesn't seem as simple as I expected.

Has anyone got any ideas how I can use sql 2008 for data storage and do away with shape files? I've read about arcsde, but really don't understand the point of another layer. No idea what it's actually supposed to do.


texts: Using IRelationalOperator ArcObjects interface with feature classes stored in Oracle?

texts:
Could you give an example how you’d use the IRelationalOperator ArcObjects interface with features from different feature classes stored in oracle database?

In my visual studio project I have an xml file which reference all my actual data and I have created variables in C# to hold the values of these data referenced in the xml file.

ESRI's example only works if you know the actual features and their index # e.g.

// Get the marsh and highway features from their respective classes.
IFeature marshFeature = vegFeatureClass.GetFeature(518);
IFeature highwayFeature = roadsFeatureClass.GetFeature(39);

// Get the geometries of the two features.
IGeometry marshGeometry = marshFeature.Shape;
IGeometry highwayGeometry = highwayFeature.Shape;

// Cast the highway's geometry to IRelationalOperator and determine if
// it crosses the marsh's geometry.
IRelationalOperator relationalOperator = (IRelationalOperator)highwayGeometry;
Boolean crosses = relationalOperator.Crosses(marshGeometry);
Console.WriteLine("Highway crosses marsh: {0}", crosses);

I want to be able to use this interface without knowing what the exact index  value of GetFeature() function is.


texts: QGIS: permanent layer of vertices is invisible

texts:
I'm working in QGIS, and I've extracted the vertices from a .shp file of polygons. If I make the vertices into a temporary layer, the points are visible and I can remove them as I please. However, this is risky because exiting QGIS means that the temporary layer gets removed.
I then tried making a permanent layer out of the vertices, but the points will not show up on my screen. Any ideas why?


texts: Map point with coordinates to zones in QGIS

texts:
I have a project in QGIS with two layers, one with several geographical zones (imported from a Shapefile file), and one with several points (imported from a CSV with longitude/latitude coordinates).

I would like to map each point to the ID of the zone that contains that point, i.e., create a new CSV in the form "longitude,latitude,ID". How can I do it in QGIS (or, possibly, any other "scripted" solution)?


texts: Fixed Buffer disappears after closing QGIS

texts:
Has anyone else experienced this problem and knows how to solve it?
When I reopen QGIS the buffer layer is still there but the buffer lines themselves are not shown in the map anymore.


texts: Convert a shapefile with northing/easting to ASCII raster with latitude/longitude

texts:
Linked below is a polygon shapefile whose projection is NAD83 and whose values are LU/LC codes. Coordinates are meters of northing and easting. I would like to convert it to an ASCII raster with coordinates measured in latitude and longitude. I'm sure it is a simple process, but I can't figure out the re-projection. Tools I have are ArcGIS or SAGA or R. I can even open the ASCII text and modify the header if no recalculation is necessary.

http://hawaii.gov/dbedt/gis/lulc.htm

EDIT: Thanks to Dan Patterson for pointing this out: the projection is actually UTM 4, and NAD83 is the DATUM. The output raster should have the same datum. If I want output to have lat/long coordinates, does that mean it has no projection?


texts: Aerial Triangulation by using GPS height

texts: Is it possible to use GPS geographic coordinates & ellipsoidal height not orthometric heights in Aerial Triangulation?


texts:
I am needing to number a selection of points following a street or set of streets using an alpha-numeric string from A001-Z999. And while I'm trying to figure out how to get these points numbered using only odd numbers (ie A001, A003, A005); more importantly, I'm needing to communicate to ArcGIS that I want them numbered starting at a certain point on the street and following the trail of selections.

I'll be needing to duplicate this process with each new, unique strand of selected points and assume each strand will travel in multiple directions. Is there a way to do this specifically within Field Calculator?


texts: Creating empty space when merging multiple GEOTIFF files?

texts:
I have multiple GEOTIFF files containing information regarding atmospheric parameters. Due to a lack of data in certain regions, I've had to do my geospatial processing separately on them.

Now I'd like to merge the files back together. When I use GDAL_merge, it creates an interpolated space between the two regions. I'd like this to remain as empty space. In addition to this, the file grows very large. Two relatively small files become many magnitudes larger.

$ gdal_merge.py -o output.tif input_1.tif input_2.tif


texts: Is it possible to use two layer.bindPopup functions in Leaflet

texts:
I have two different layers in a leaflet map. I have one layer in which I would like the popup to return some text from the layer's attributes. They other layer I would like the popup to return an image which I have embedded the link to in the layers attributes. I have successfully implemented both of the layer.bindpopup functions successfully individually but I'm having no luck implementing them together.

function onPlot(feature, layer) {
    if (feature.properties) {
        layer.bindPopup(
            '');
    }
};

function onEachFeature(feature, layer) {
    if (feature.properties) {
        layer.bindPopup(
            '' + feature.properties.Anum +
            '
 Hunting District: ' +
            feature.properties.DISTRICT +
            "
 Elk Hunter Success " +
            feature.properties.Success +" %");
    }
};

function style (feature,layer){
    if (feature.properties) {
        return  {
            color: feature.properties.color,
            fillColor: feature.properties.fillColor,
            fillOpacity: feature.properties.fillOpacity
        };
    }
};

var Plot= new L.geoJson(Plots, {
    onPlot:onPlot
    style:style
}).addTo(map);

var S_2006= new L.geoJson(Success_2006, {
    onEachFeature: onEachFeature,
    style:style
}).addTo(map);


texts:
I have point feature with name PointPhoto.
This feature have attribute:

id
area
PhotoPath

Now I am using GeoExt popup to display this attributes, but I can't find a way to display image from image path, this is the code:

 <html>
    <head>
        <link rel="stylesheet" type="text/css" href="ext/resources/css/ext-all.css">
        <script type="text/javascript" src="ext/adapter/ext/ext-base.js"></script>
            <script type="text/javascript" src="ext/ext-all-debug.js"></script>
        <script src="../OpenLayers-2.9.1/OpenLayers.js"></script>
        <script type="text/javascript" src="geoext/script/GeoExt.js"></script>
          <link rel="stylesheet" type="text/css" href="http://localhost:8181/geoserver/openlayers/theme/default/style.css"/>
        <script type="text/javascript" src="../proj4js/lib/proj4js-compressed.js"></script>
        <link rel="stylesheet" type="text/css" href="geoext/resources/css/popup.css">
        <script type="text/javascript">
var app, items = [], controls = [];
var store;
Ext.onReady(function() {
            app = new Ext.Viewport({
                layout: "border",
                  items: items
            });
        });
        items.push({
            xtype: "gx_mappanel",
            ref: "mapPanel",
            region: "center",
            map: {
                numZoomLevels: 19,
                controls: controls
            },
            extent: OpenLayers.Bounds.fromArray([
                    2189574.7042680467, 5418655.108865695,
                    2199298.319819289, 5426457.253025988
                    ]),
            layers: [new OpenLayers.Layer.WMS(
                "cite:PointPhoto",
                "http://192.168.1.2:8181/geoserver/cite/wms",
                {layers: "cite:PointPhoto"},
                {isBaseLayer: false}
            )] ,
            });
            controls.push(
        new OpenLayers.Control.PanZoomBar({position: new OpenLayers.Pixel(2, 15)}),
                new OpenLayers.Control.Navigation()
           );

controls.push(new OpenLayers.Control.WMSGetFeatureInfo({
autoActivate: true,
infoFormat: "application/vnd.ogc.gml",
maxFeatures: 3,
eventListeners: {
        "getfeatureinfo": function(e) {
            var items = [];
            Ext.each(e.features, function(feature) {
                items.push({
                    xtype: "propertygrid",
                    title: feature.fid,
                    source: feature.attributes
                         });
            });

            new GeoExt.Popup({
                title: "Feature Info",
                width: 200,
                height: 200,
                layout: "accordion",
                map: app.mapPanel,
                location: e.xy,
                items: items[0]
            }).show();
        }
    }
}));

     </script>
    </head>
    <body>

    </body>
</html>

My question is how I can display image inside of this popup?


texts: QGIS not asking to save my modified MapInfo TABfile

texts:
Editing TAB files works fine from user interface and from python code too. If user click SAVE button from QGIS or from python code - it works FINE. BUT if user is closing directly QGIS application - there is not question to save changes in TAB file :( For shape files QGIS is asking before application shut down.

Any ideas?


texts:
I've install GeoNode through apt-get install geonode. As per the Quick Installation Guide.
After such installation how do one customize the pages ? I don't seem to find any indication reading the doc on how it should be done.

I followed the instructions at Setup & Configure HTTPD and I'm having some issues customizing.
The changes made in site_index.html don't show up when I visit the home page.

Am i obligated to create a project as per GeoNode (v2.6) installation on Ubuntu ?

Any advice would be nice.


texts: How to visualize data for a specific class in change detection plot in R?

texts:
I have created a change detection plot shown in image below for two rasters having 6 classes each using the code below. Now the plot, shows the change for all the classes but i am interested in only visualizing the change for class 4 and 5 i.e. pixels where class 4/5 changed to another class or where another class changed to 4/5.

How can i do this task?

Here is my sample code:

library(raster)
library(rasterVis)

r <- raster()

set.seed(123)

lc1 <- setValues(r, sample(1:6, 64800, replace = T))
lc2 <- setValues(r, sample(1:6, 64800, replace = T))
lc1_uniq <- unique(lc1)
lc2_uniq <- unique(lc2)
grid_ <- expand.grid(lc1_uniq,lc2_uniq)
names(grid_) <- c('from','to')
grid_$code <- 1:dim(grid_)[1]
grid_$change <- grid_[,1] != grid_[,2]

head(grid_)
change <- function(x){
  grid_[x[1] == grid_[,1] & x[2] == grid_[,2],'code']
}
changeDet1 <- calc(stack(lc1,lc2), fun = change)
codes_ <- data.frame(ID = grid_$code,value = paste0('from ',grid_[,1],' to ',grid_[,2]))
logical_test <- which(grid_$change == T) # remove no change classes
codes_ <- codes_[logical_test,]
# Create a Raster Attribute Table
rat <- levels(changeDet1)[[1]]
rat[["Changes"]] <- codes_
rat[["Changes"]] <- c("land","ocean/lake", "rivers","water bodies")
levels(changeDet1) <- rat
# Plot
levelplot(changeDet1, par.settings=PuOrTheme(), xlab="", ylab="")

This question continues Converting output of crosstab() in R to raster; from classified raster images to show overall land cover change in different classes?. All the details (plots,comments) might help in understanding this question more.


texts: Buffer tool not surrounding complete shape in ArcMap

texts:
I'm running a basic buffer tool on a polyline shapefile (teal line), to create a 50km buffer (beige) around it. However, ArcMap appears to only buffer part of the map (see image below).

How can this be worked around? Is the polyline too complicated?

I have already used the simplify line tool to try and speed up the process.

Storage type of buffer target: Shapefile

Co-ordinate system: WGS84

Co-ordinate resolution: Initially 1:250,000, but run under a 'Simplify Line' to a simplification tolerance of 0.5km

Desktop Version: 10.2


texts: Calculating Volume of Reservoir using python and ArcMap10?

texts: I have a raster file with different water depth in each cell. I am using Aster DEM so my resolution of that file is 30m by 30m. Now I want to calculate the total volume in that area by multiplying water depth in each cell by area of the cell and calculating the total. How can i perform this using python and arcmap 10 ? I would be glad to receive your answer


texts:
i have a problem with annotation (ok i will explain) :
in database i have :

parcelle : a geometry table that contain polygons
personne : a table contain parcel's owners
pll_personne : is the relationship (n,n) between the (1) & (2)

what i want is : to display parcels with annotation Name of there owners see picture please

there is somebody that can help i'm very Jam !! :(


texts: Clipping polygon to parcel without affecting polygons in neighboring parcels using ArcMap?

texts:
In one of my weekly QA tasks using ArcMap 10.3, I receive a table of parcel changes that have occurred. I need to review the parcels for polygons that overlap the parcel line. I use a tedious combination of reshape, cut, and clip editing tools. Using Clip is the quickest but here's my problem:

The polygons are snapped to the parcel line. When I clip a polygon to the parcel with "Preserve the area that intersects" selected, another polygon in the neighboring parcel that shares the same boundary will then be deleted.

I want to know if

a) Outside of code, is there a quicker way to do this without deleting neighboring features? Perhaps within ModelBuilder?

b) How simple would it be to write a code to bulk clip selected features based on their corresponding ParcelID ?


texts: I need to plot the field lines of a 3D vector (Vx, Vy, Vz) of which I have the 3 components at each point of a given 3D grid. One way is to use FLOW3. However, in order to use FLOW3 I need to have the vector in a form of 3D array, and I don't have it. Instead I have a table of six columns in the form X,Y,Z,Vx,Vy,Vz: at each row I have the spatial position (X,Y,Z) and the three components of the vector (Vx,Vy,Vz) at that position. Can REFORM or MAKE_ARRAY do what I need? I can't find a simple way to do it.


texts: How to use the multiprocessing module on python to execute a spatial join?

texts:
I have several shapefiles of the entire Amazon rainforest each with millions of polygons (they are fairly detailed and take some time to draw on ArcMap). I need each of these shapefiles to be spatially joined with another shapefile called landuseA (which, itself, is a very big file). Given how computer intensive this task is, I would like to use the multiprocessing module on python to improve performance. How do I do this?

I have read the blog post by Esri on the multiprocessing module, but I cannot figure out how I can use their example to write a script that will use multiple processors to execute a spatial join. Blog post is here: http://blogs.esri.com/esri/arcgis/2012/09/26/distributed-processing-with-arcgis-part-1/

If you have other suggestions for how this task might be better accomplished, feel free to share that too.


texts: I am trying to display a map/shapefile of an NTS250k grid of Canada as shown in MS Power BI. BI requires to import the shapefile using .json format. I am just wondering if there is a way to be able to overlay a map of Canada onto this (without affecting the shape selection) so the user can see what it is they are looking at.


texts:
I have a CAD drawing that opens in true view (the free version) and I identified the layer I wish to export. Now, I want to bring the layer into ArcGIS and merge it into my geodatabase. However, the layer does not show up in ArcMap. In fact, I only get one layer in arcmap with actual objects in the table or features drawn on the map: xref. aka the base of the drawing.

So far I have tried turning all drawing layers on and off, adding a spatial reference, and converting CAD to GDB using the gp tool in ArcMap. Here is a screenshot of the layer in TrueView:

Here is a screenshot of what is in ArcMap:

My question: How can I get this layer into ArcMap with only the free CAD resources available to me? How can I even diagnose the problem? I am unsure where to begin..


texts:
I installed PostgreSQL 9.2.4 and PostGIS 2.0.3 on my Windows 8 machine. Then I created a table using phpPgAdmin and added a geometry column using the following statement.

SELECT AddGeometryColumn('tablename', 'area', -1, 'POLYGON', 2);

That resulted in these table columns.

+------+---------------------+--------+----------------------------------+-----------+
|Column|Type                 |Not Null|Default                           |Constraints|
+------+---------------------+--------+----------------------------------+-----------+
|id    |integer              |NOT NULL|nextval('simple_id_seq'::regclass)|PRIMARY KEY|
|name  |character varying(50)|        |                                  |           |
|area  |geometry(Polygon)    |        |                                  |           |
+------+---------------------+--------+----------------------------------+-----------+

After some simple text queries I tried to play with geometry data. This is the query I tried to execute.

INSERT INTO tablename (name, area)
VALUES (
    'Lake',
    GeometryFromText('POLYGON((84 104, 212 48, 321 122, 243 179, 275 238, 131 240))', -1)
);

Unfortunately, it throws an error message.

ERROR:  function geometryfromtext(unknown, integer) does not exist
LINE 4:     GeometryFromText('POLYGON((84 104, 212 48, 321 122, 243 ...
            ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Since AddGeometryColumn was successful, I assume that PostGIS is set up properly. How can I insert geometry data in my table? What am I doing wrong?


texts: Getting Started Using Leaflet Draw to CrowdSource Shapes on a Map

texts: I want to crowdsource shapes on a map for an organization I'm working with. It seems like Leaflet Draw something along the lines of the Mapping Neighborhoods projects. As someone with some experience both working with spatial databases, and making Leaflet maps, I was wondering if there's a tutorial available to get started combining a Leaflet map writing to a database.


texts: ArcGIS - after merging multiple shapefiles, there is a boundary between them I want rid of. How do I get rid of them?

texts:
I have downloaded four geology tiles and merged them together. It works, and I can then go on to customise the tiles to colour code each geological unit no problem, but those four lines joining the tiles I merged look so messy and I want rid of them. I get these lines when I merge shapefiles to make the merged file, and when I merge layers to make the merged file. No solution as of yet.

I imagine it could be a stylistic thing, where I need to make the borders transparent or something. Thoughts?


texts: comparison of Historical aerial photographs to today's satellite images

texts:
I'm working on a project in which I have to compare vegetation cover changes. So we have 1950's aerial photographs which must be compared to actual satellite images. Images must be properly adapted for the satellite image, so I can measure the "before" and "after" accordingly. Do you all know any tutorial that teaches how to do this? I'll use Quantum GIS for these procedures.

The area I'm working is a landscape of ca. 50.000 ha. The aerial photographs were already digitized and yes, they're black-white. So I have several images, each corresponding to a certain amount of landscape. These photographs are not georeferenced, so I'll have to locate their coordinates manually. This is not a hard task, the hardest part for me is to properly deal with the errors / or distortions, photography angle, etc. in order to match to the satellite images. I think the resolution is quite OK to work with, since my idea is to verify how much vegetation advanced into grasslands in the past 60 years.


texts: Using Field Calculator in QGIS to update columns with centroid coordinates

texts:
In QGIS I would like to update two columns with the centroid easting and northing of my polygons. I am using the field calculator. I have tried using $x and $y but this only seems to work for points.

I have checked the help but it does not seem to mention being able to do this.


texts: How to update kml file/layer dynamically?

texts:
I am using a kml file for displaying map with custom boundary lines. That kml file is downloaded from some website. In that file, inside the placemark tag there is no point tag for display icons. Eg:

<Placemark>
    <name>Spot 2</name>
    <description>.....</description>
    <styleUrl>....</styleUrl>
    <MultiGeometry><Polygon><outerBoundaryIs><LinearRing><coordinates>
    .........
    </coordinates></LinearRing></outerBoundaryIs></Polygon></MultiGeometry>
</Placemark>

That is the placemark tag contains. 1) How can I add point tag for all placemark tag. Is there any way to add dynamically? My kml file has 5000 and above placemarks. 2) Point tag coordinates will refers the center of the polygon.

i.e) I need the following

<Placemark>
    <name>Spot 2</name>
    <description>.....</description>
    <styleUrl>....</styleUrl>
    <MultiGeometry><Polygon><outerBoundaryIs><LinearRing><coordinates>
    .........
    </coordinates></LinearRing></outerBoundaryIs></Polygon></MultiGeometry>
    <Point>
        <coordinates>144.253,-36.6632,0</coordinates>
    </Point>
</Placemark>

Note:

Iam using geoxml3 parser to display kml layer in google map.


texts:
I have a dataset (vector data) that contains the polygons and their corresponding crop type. This data is highly imbalanced i.e, I have a lot of classes that contain one sample only, as shown in the photo down below. To oversee this issue I tried to implement SMOTE but this throws an error because SMOTE needs at least 2 samples of the same class to create new ones (k_neighbors=1 by default).
I can't just use the dataset as it is because the results are going to be biased towards large classes and (2) some algorithms throw an error during training.

What I don't understand is that the minority class (correspond to one polygon in the vector data) would cover more than one pixel (when matched with the features, I have 15 features) therefore, at least 15 pixels correspond to one of the minority classes.
This the error I get:
ValueError: Expected n_neighbors <= n_samples,  but n_samples = 1, n_neighbors = 2

Is there a way to overcome this issue? Should I duplicate the existing samples and then use SMOTE to oversample them?


texts: Error while performing line segment detection in monteverdi

texts:
I am not able to perform Line Segment Detection in this stand alone tool 'Monteverdi'. I have already tried it number of times but no output generated and got the error
''LineSegmentDetection' failed with return status -1', please help or suggest some other way to perform this task.

result og log panel plugin tab :

2017 Oct 28 12:53:56 : Application.logger (INFO)
otbcli_LineSegmentDetection.bat -in F:/Results_Monteverdi/GE8_Smooth.tif -out F:/Results_Monteverdi/LSD_ge8 -ram
128 2017 Oct 28 12:53:58 : Application.logger (FATAL) The following error occurred during OTB-application execution:
c:\dashboard\otb\src\modules\io\vectordataio\include\otbVectorDataFileWriter.txx:138: Could not create IO object for file F:/Results_Monteverdi/LSD_ge8 Tried to create one of the following: OGRVectorDataIO You probably failed to set a file suffix, or set the suffix to an unsupported type

I was trying the Line Segment Detection for GeoTIFF images and after trying a lot was not able to get the output. Later i converted the input images to PNG format and removed the properties and personal information contained in the image, after that only LSD was able to process but now the output which is generated by the tool , is like the mirror image of the line segments. I am trying to rectify but if anyone who already tried this is welcome for the help. I have attached the images below:


texts: Splitting lines at a specified distance and add attribute

texts:
I have used QChainage to create points along a line every meter. What I am ultimately looking for, is to split the line layer along this point into segments with the chainage distance as an attribute. I use QGIS 2.6, but am open to options using GRASS standalone. I have installed the LRS plugin, but would like to know more about how to use to take me closer to my goals.

The QChainage gives me points on the line where these numbers are :

0-----1-----2-----3-----4-----5

I would like to split the line into segments where the QChainage points are and add the distance attribute to the segment.

|--1--|--2--|--3--|--4--|--5--|

I started tinkering with the GRASS v.lrs functions, but have not had success with it. I'm blocked by an error:

WARNING: Not enough points attached to the line


texts: Search for string and only display characters before it in Mapinfo label expressions

texts:
So I have a property table that has addresses like below.

12 Andrews CL, Pxxx
1/12 Andrews CL, Pxxx
2/31 Macrossan ST, Pxxx
2/7-11 Barrier   ST, Pxxx

I need to only show

12
1/12
2/31
2/7-11

Can this be done within the expression builder in Mapinfo 12?

In python it would be something like address[:address.find(' ')]


texts: Creating individual offline layers for ArcGIS Runtime (100.x) map

texts:
I'm wanting to convert a fairly large ArcGIS Runtime 10.2.7 WPF project to 100.x   The project uses entirely offline content, which is stored in Runtime Geodatabases that are created in ArcMap (File ==> Share As ==> ArcGIS Runtime Content)  I programmatically build the map using the layers in multiple runtime geodatabases.  The map is entirely read-only, there is no data editing need.
At this website https://developers.arcgis.com/net/offline-maps-scenes-and-data/  about halfway down on the page, it says the following:
An alternative to taking a map offline is to construct the map using individual offline layers. Layers can either reference offline data sources on the device or be exported from suitable online services. Working directly with layers gives you full flexibility to compose your own offline map.
This is exactly what I want to do (using data stored on the device), but of course, it doesn't give any recommendation or code samples on how to do it.  Basically what I'm looking for is what the best way is to package the feature layers that I want to use in a way that's best for a custom (programmatically) made offline map.
What I've tried so far...In ArcGIS Pro 2.7

Package the data in a mobile map package.  I've sort of gotten this to work, but you have to "remove the layers" from the mobile map first before you can use them in your own map. Not sure if mobile map package is the ideal solution because I don't want to use the map in the mobile map package (I just want to use the feature (data) layers)

Create a mobile geodatabase (under Data Management Tools ==> Create Mobile Geodatabase) and then add my feature layers to that (using the "Feature Class to Geodatabase" tool)  However, when I open the mobile geodatabase in a WPF app, the Geodatabase.GeodatabaseFeatureTables list is empty (it's like it doesn't recognize that my feature layer as a GeodatabaseFeatureTable)

I need to be able to export/package my feature layers using the "ArcGIS Pro version of ArcPy" (preferably this instead of using the old Python 2.7 version...)  So this eliminates the "Create Runtime Content" method used in ArcMap.
Thoughts/ideas/suggestions welcome.


texts: Proccessing tool to add a new field to existing vector attribute table

texts:
Isnt there a tool to add a new field to a existing vector table without creating a new dataset, like Add Field in ArcGIS? All i can find is "Add field to attributes table" which creates a new vector dataset:

  Add field to attributes table This algorithm adds a new attribute to a
  vector layer. The name and characteristics of the attribute are
  defined as parameters. The new attribute is not added to the input
  layer but a new layer is generated instead.

I know it is possible in the attribute table by editing - New field, but I want to create a model which starts with an empty Point layer and adds 15 fields to it.


texts: Trace Geometric Network in custom script tool?

texts:
I'm trying to make a Python script tool which uses the Trace Geometric Network tool (10.3). My script works in the Python window, but as a script-tool, it runs and NOTHING HAPPENS. Messages state: Executing..., Start time..., Running script..., Completed script... Succeeded at... It acts like the tool executed with no errors, but there's no visible output.

The Trace Geometric Network tool refuses to work inside a Python script tool the way it does in the Python window.  I tried playing with the parameters, data types in of parameter settings, and deleting and remaking the tool several times. I even tried using other geometric networks with the tool. Finally, I made a bare-bones script with only the trace tool to troubleshoot, and NOTHING. I get the same action as above.

Here's my bare-bones script:

geoNetwork = arcpy.GetParameterAsText(0)

inflags = arcpy.GetParameterAsText(1)

barriers = arcpy.GetParameterAsText(2)

arcpy.TraceGeometricNetwork_management(geoNetwork, "Trace", inflags, "TRACE_DOWNSTREAM", barriers)

Here's a screen shot of how I'm setting my parameters in the tool properties:

How can I to get this to work as a Python script tool?

Here's the actual code. As this tool executes, I can see it looping over the flags and selecting them one by one. No output group layer files are added to the map and no apparent traces are made. Again, these functions work as intended if pasted in the python window, but misbehave as python-script tools.

inflags = arcpy.GetParameterAsText(0)
geoNetwork = arcpy.GetParameterAsText(1)
barriers = arcpy.GetParameterAsText(2)
outPath = arcpy.GetParameterAsText(3)
outName = arcpy.GetParameterAsText(4)

def makeFullPath(path, name):
    return path + "\\" + name

def TraceCapture(inflags, geoNetwork, barriers, outPath, outName):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    flagList = list(arcpy.da.SearchCursor(inflags, 'UNITID'))
    trace = 1
    segmentCapture = []
    for flag in range(1, len(flagList)+1):
        try:
            outputLayer = "Trace_" + str(trace)
            expression = '"OBJECTID" =' + str(trace)
            arcpy.SelectLayerByAttribute_management(inflags, "NEW_SELECTION", expression)
            arcpy.TraceGeometricNetwork_management(geoNetwork, outputLayer, inflags, "TRACE_DOWNSTREAM", barriers)
            arcpy.AddMessage("Trace_" + str(trace) + " complete")

            for lyr in arcpy.mapping.ListLayers(mxd, "", df):
                if lyr.name == "Sewers":
                    pipesCursor = arcpy.da.SearchCursor(lyr, 'ASSETID')
                    for i in pipesCursor:
                        segmentCapture.append((flagList[trace - 1][0], i[0]))
                if lyr.name == outputLayer:
                    arcpy.mapping.RemoveLayer(df, lyr)
            trace += 1
        except RuntimeError:
            print "Flag " + str(trace) + " was invalid and skipped."
            trace += 1
            continue

    arcpy.CreateTable_management(outPath, outName)
    outFile = makeFullPath(outPath, outName)
    arcpy.AddField_management(outFile, "manholeID", "TEXT", "", "", 10)
    arcpy.AddField_management(outFile, "pipeID", "TEXT", "", "", 10)

    with arcpy.da.InsertCursor(outFile, ["manholeID", "pipeID"]) as cursor:
        for row in segmentCapture:
            cursor.insertRow(row)

TraceCapture(inflags, geoNetowrk, barriers, outPath, outName)


texts: swipe dynamic layers using spy glass in ol3

texts:
I am working on swipe based ol3 web application where i want to use spy glass for dynamic added layers.
I have seen layer spy example but this is only for single layer swipe ..
I want for dynamic layers.


texts: Why are some columns not recognized as numeric even though I have a CSVT?

texts:
I have problems importing a .csv file.

Some Values of the file are not supposed to be formated in text style, but in number style, thus i created an additional .csvt file like explained here.

Both files with exactly the same name and the same location are created. The number of columns is the same and I add the .csv file with the "Add Delimited Text Layer" tool. Nevertheless I am not able to choose any columns for a graduated symbology.

My .csv file includes 12 rows, of which the 5th, 11th, and 12th column consist numbers.
My .csvt file therefore looks like that:

"String","String","String","String","Integer","String","String","String","String","String","Integer","Integer"

I tried to create it with the standard editor from Windows and Notepad++.
Are there any common mistakes done by newbies, like I am one?
Thank you for your help!


texts: Parallel lines not showing on my map document (QGIS; v.parallel function)

texts:
Weird issue: I am trying to use v.parallel to draw a line parallel to a single line layer I have. It processes without issue, but the parallel line does not appear on the map. The curious thing is that I can click zoom to layer and it zooms to the area where the line SHOULD appear...it's just not physically there.
These are the parameters I set. Why is it not showing?

Also it has no apparent ability to loop through several lines if I wanted say, 10 parallel lines at an equal spacing. Do I really have to manually do each parallel line?


texts:
I do not know what  happened but I have lost all my ArcToolbox log history that was available trough "Results" menu.

The strange thing is that I can access the log files that are saved on my harddrive. Those files are on .xml format but I do not know how to import them back to the "Results" history menu. In this picture you can see that I have very few logs in the "Results" history, but I have all of the history files!

Is there any way I can import back these files to the "Results" menu?

Can someone give me a help on this? To me is extremely useful to see the results history! It helps me a lot!


texts: Converting point shapefile with N points to N layers with one point each using Global Mapper

texts:
I have a point shapefile layer that contains 50 points. I want to create 50 layers with these points. Meaning every point ends up in its own layer.

How can I do this using Global Mapper?


texts: Error loading GeoJSON data - XMLHttpRequest cannot load file?

texts:
When I try to load a geojson format vector data, I get the following error :

XMLHttpRequest cannot load file. Cross origin requests are only supported for protocol schemes: http, data, chrome, chrome-extension, https.

But when I start a simple python server using http.server, everything seem to work. What's the problem?


texts: Data is not appearing in a local OSM server after applying changes to it

texts:
Our company has a local OpenStreetMap server hosted on its own server and domain. Our cartographers draw streets and houses there in edit mode.

The problem is that the houses are visible only in the edit mode, and in the view mode they do not appear at all. I know that in normal OSM it takes several hours before the changes are published. But on this local server, the changes have not been published for two weeks now.

Here is an example. There are no houses and streets in the view mode:

In the local server's edit mode everything is there:

(Maybe it's important, but in the normal OSM's edit mode there are no houses and streets, too).

The employee who created the local server quit and does not get in touch, I'm new to OSM. Please advise me where to look for a solution, in which direction?


texts: Measuring distances between multiple polygons whilst avoiding buildings?

texts:
I am trying to map habitat corridors for different species of animals.

A map with both habitat areas and buildings as shape files (shown in the image below)
Maximum distance an animal will leave its habitat for
Maximum distance an animal will travel

I want to map which habitat areas will link up.

I tried Cost Distance but I am not factoring terrain change in elevation so I don't think this will work.
I have also done euclidean distance with the habitat areas but I can't figure out how to get it to recognize buildings - at the moment it just goes through them.


texts: Cluster Tooltip to display list of features

texts: Trying to work out how to make the tool tip show a list of features represented by a Point Cluster. When I use identify on a Point Cluster it provides a list of the features in that cluster to identify. Just trying to make the same list appear on mouse over. QGIS 3.4


texts: What's the usual term for "ground plan based"

texts:
I don't know if "ground plan based" is the right word in the following context:

A computer program reads a ground plan and produces a result ground plan. The result groundplan's deviation from the input ground plan must be within a prescribed tolerance to be called "ground plan based".

Possible words for this situation might be:

ground plan based, topographic planimetry, planimetric reconstruction

But I'm not sure which is the common word? Which word would you use to describe that both of the ground plans match nearly 100%?


texts: NameError from using definition query in ArcPy?

texts:
I am having trouble getting my definition query to be defined. What am I missing? The current error is:

Runtime error
Traceback (most recent call last):
  File "<string>", line 44, in <module>
NameError: name 'df' is not defined

# Import arcpy module
import arcpy
import os

# Local variables:
plain_radius5Mi_lyr = "Filepath\\plain_radius5Mi.lyr"
hike__2_ = plain_radius5Mi_lyr
XYplain2__2_ = "XYplain2"
XYplain3 = "Filepath\\plain.gdb\\XYplain3"
arcpy.env.overwriteOutput = True
hike = "Filepath\\plain.gdb\\hike"

# Process: Select
arcpy.Select_analysis(XYplain2__2_, XYplain3, "")

# Process: Multiple Ring Buffer
arcpy.MultipleRingBuffer_analysis(XYplain3, hike, "5;10;15;20;25;35;50", "NauticalMiles", "distance", "NONE", "FULL")

# Process: Apply Symbology From Layer
arcpy.ApplySymbologyFromLayer_management("hike", plain_radius5Mi_lyr)
# top level directory where all the map docs are
path = r"<path>"

# recursively search the directory for .mxd files
def getMXD(path):
    for dirpath, dirnames, filenames in os.walk(path):
        for f in filenames:
            if f.endswith(".mxd"):
                yield os.path.join(dirpath, f)

for mapdoc in getMXD(path):
    # might need to change wildcarding here
    mxd = arcpy.mapping.MapDocument(mapdoc)
    df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
for lyr in arcpy.mapping.ListLayers(hike, "*", df):
     if lyr.supports("DEFINITIONQUERY"):
          lyr.definitionQuery = "ID = '35'"


texts: How to calculate angle between two points and vertical

texts:
Is there an easy method to calculate angle between two points and vertical with PyQGIS ?
In QGIS I use Cad-Tools extension, I select two points and it give me an azimuth. Now I would like to do the same in my script where i'm already able to select two points.


texts: Viewshed analysis plugin fails with 'ascii' codec can't encode characters

texts:
I have problem with the Viewshed analysis plugin for QGIS. After I set all the options in it and I run it, every time I get this messeage:

An error has occurred while executing Python code:

UnicodeEncodeError: 'ascii' codec can't encode characters in position 0-3: ordinal not in range(128)
Traceback (most recent call last):
  File "C:/Users/Dzole/.qgis2/python/plugins\ViewshedAnalysis\viewshedanalysis.py", line 93, in load_cmbObsField
def load_cmbObsField(self,cmb_index): self.reload_dependent_combos(cmb_index,'cmbObsField')
  File "C:/Users/Dzole/.qgis2/python/plugins\ViewshedAnalysis\viewshedanalysis.py", line 109, in reload_dependent_combos
ly_name=str(l)
UnicodeEncodeError: 'ascii' codec can't encode characters in position 0-3: ordinal not in range(128)

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.14.0-Essen Essen, 98e9610

Anyone knows what's this about?


texts: What do these curly graphics represent on this map?

texts:
What do the two curly graphics to the left and right of this map represent? They are the only two on the map and seem to be attached to the border of an area.

I'm sorry I haven't posted the full map, I'm not sure I'm allowed.


texts: Find closest point between polygons inside other polygon

texts:
I have access to ArcInfo (ArcGIS 10.1). The dark polygons are settlements. The green polygon is derived from a suitability analysis. The surrounding white area is restricted.
I want to find a point that will serve as a bus stop that has the shortest possible distance between the 2 settlements but it must fall within the green area.
The initial idea was to convert the settlements to raster, raster to point (to mimic density), merge points and mean center. In this case we are not taking into account the green area, so there is a great possibility that the mean center will fall inside the bigger settlement. However I want it to be specifically inside the green polygon. Do you have any ideas?


texts: Get DISTINCT nearest neighbors with indexed search

texts:
I'm using the indexed nearest neighbor search to return the closest points to a given starting location. My points also have attributes like color and shape, stored as strings. I want each query to only return the closest point for each type of point.

So if my points were, in order of distance to my starting location:

red square
red square
red circle
orange square
blue circle

The query for the four closest points would return:

the closest red square
red circle
orange square
blue circle

I tried using DISTINCT but it's asking to give the DISTINCT column in the ORDER BY and  I haven't had success in doing that.


texts: Is it possible to find which sectors of a site are pointing to a specific lat long at a distance away?


texts:
I have a two tables.
1 table have srid=4326
2 table have srid=70066

I want to compare geometry columns of this tables. I find function ST_Equals(geometry, geometry) but its not what i want. I need something what can compare geometries which slightly different from each other like on picture.

UPDATE
i have test polygon in tables with area=9098.79m^2 in first table and with area=9721.82m^2 in sacond table.
I try this query

SELECT   st_area(ST_Intersection(ST_TRANSFORM(filedata.the_geom,70066),filedata_temp.the_geom))
FROM filedata, filedata_temp
Where filedata.num=1

srid=70066 is rectangular coordinate system

and i get query result 8910.5625
why so much?


texts: Is there a faster way to add a field based on zonal statistics to an attribute table?

texts:
Let's say I have a shapefile 'Mypolygons' with an ID field and a raster 'DEM' giving elevation values. I want to add a field for each row in the 'Mypolygons' attribute table to give the average elevation within each polygon.

My current approach is the following:

Add field "Avg_Elev" to 'Mypolygons'.
Use Zonal Statistics as Table with 'Mypolygons' as the zone data, 'ID' as the Zone field and 'DEM' as the in_value_raster, and 'Output' as the out_table to hold mean elevation values.
Join 'Mypolygons' and 'Output' based on the 'ID' field.
Use Calculate field to populate values of 'Mypolygons' 'Avg_Elev' field based on the 'MEAN' field in 'Output'.
Remove join.

Is there a single tool or a more efficient way to perform this task? A python solution will work for me as well.


texts:
Trying to convert a CSV file (3 columns: Lon, Lat, Ref) to .tif GeoTIFF (the Ref values) using gdal. Using Python 3.6, I am getting the error below.

What am I doing wrong?

C:\ProgramData\Anaconda2\Envs\Python3.6\lib\site-packages\osgeo\gdal.py in Grid(destName, srcDS, **kwargs) 966 srcDS= OpenEX(srcDS, OF_VECTOR) 967 968 return GridInternal(destName, srcDS, opts, callback, callback_data) 969 970 def RasterizeOptions(options = [], format = None,

c:\ProgramData\Anaconda2\envs\Python3.6\lib\site-packages\osgeo\gdal.py
in GridINternal(*args) 3227 def GridInternal(*args): 3228 """GridInternal(char const * dest, Dataset dataset, GDALGridOptions options, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset """ 3229 return_gdal.GridInternal(*args) 32230 class GDALRasterizeOptions(_object): 3231 """Proxy of C++ GDALRasterizeOptions class.""" ValueError: Received a NULL pointer.

Code:

import os
import gdal
#this code will loop over a directory and read all .csv files
dir_with_csvs = r"C:\User"
os.chdir(dir_with_csvs)

def find_csv_filenames(path_to_dir, suffix=".csv"):
    filenames = os.listdir(path_to_dir)
    return [ filename for filename in filenames if filename.endswith(suffix) ]
csvfiles = find_csv_filenames(dir_with_csvs)
for fn in csvfiles:
    vrt_fn = fn.replace(".csv", ".vrt")
    lyr_name = fn.replace('.csv', '')
    out_tif = fn.replace('.csv', '.tiff')
#create a .vrt file
    with open(vrt_fn, 'w') as fn_vrt:
        fn_vrt.write('<OGRVRTDataSource>\n')
        fn_vrt.write('\t<OGRVRTLayer name="%s">\n' % lyr_name)
        fn_vrt.write('\t\t<SrcDataSource>%s</SrcDataSource>\n' % fn)
        fn_vrt.write('\t\t<GeometryType>wkbPoint</GeometryType>\n')
        fn_vrt.write('\t\t<GeometryField encoding="PointFromColumns" x="Lon" y="Lat" z="Ref"/>\n')
        fn_vrt.write('\t</OGRVRTLayer>\n')
        fn_vrt.write('</OGRVRTDataSource>\n')
#create 2 .tif files
output = gdal.Grid('outcome.tif','name.vrt')

output2 = gdal.Grid('outcome2.tif','name.vrt', algorithm='invdist:power=2.0:smoothing=1.0')

CSV sample:
Lon,Lat, Ref

-113.3869,52.8681,30.5
-113.3718,52.8771,0.5
-113.5384,52.8951,6.5
-113.6292,52.9131,0.5
-113.6141,52.9131,0.5
-113.5989,52.9311,15.0
-113.5838,52.9311,0.5
-113.7049,52.9401,0.5
-113.5535,52.9401,1.5
-113.6898,52.9491,0.5
-113.4021,52.9581,22.0
-113.5384,52.9761,0.5
-113.5535,52.9851,0.5
-113.5384,52.9851,0.5
-115.7797,53.0031,0.5

This is my geotiff. Pixels image are just an example and not at same time of mygeotiff file (don't expect a match).

so here is what i did: import subprocess (below import gdal), removed the last 2 lines of code (below #create 2 .tif files), and replaced them with (indendent):

output = 'gdal_grid -a invdist:power=2.0:smoothing=1.0 -zfield "Ref" -of GTiff -ot Float64 -l %s %s %s' % (lyr_name, vrt_fn, out_tif)
subprocess.call(output)

The output was not what I expected at all. my aim was to create a colored pixels (centroid is Lon-Lat and color scale is based on "Ref" value). Can I still do that with my geotiff file (below)?


texts: Popup that contains feature type details from geoserver

texts:
I would like to seek for help regarding my project. I'm new in Web Mapping. And as of now I have WMS layer (which is polygon type) from geoserver that is overlayed in a WMS openlayers basemap.

Now I would like to have an interaction, a popup to be specific, that displays the feature type details that is from the geoserver. Also before that the user could highlight the selected polygon or feature.

Please enlighten me on what first to prepare and what next to do.


texts: Statistical Analysis of Shapefile Proximity within Raster Image using ArcMap?

texts:
I have a raster image that has been reclassified into 15 equally incremented values/levels. This raster represents "sinkhole risk" in a certain county. So level 1 is the lowest risk region and level 15 is the highest risk region.

I also have the shapefile data for the actual sinkholes existing in that county.

My hope is to come up with some method of statistically analyzing the data by answering the following question:

"What percentage of the sinkholes are in the risk regions 1-3?"
or

"What percentage of the sinkholes are in the risk regions 4-6?"
etc...

Any suggestions on a tool to use or a general approach? I just need some guidance in the general right direction!

I posted a photo below for more clarity! The photo has been simplified into 3 values/levels, but my real map has 15 values.

I am using ArcMap 10.1.


texts: How to plot a point on the highest elevation of a raster on a polygon. QGIS

texts:
I have a shp with more than 1700 polygons on a elevation raster.

I need to plot a point that represents the highest altitude location within the polygons, and then export those points as a new shp.


texts: exec(open('/tmp/tmpu7edmqsu.py'.encode('utf-8')).read()) on PyQGIS 3

texts:
I tried to load csv layer using this script:

file = "/home/ayato/Desktop/26mei.csv?delimiter=%s&xField=%s&yField=%s"
uri = file % (",", "x", "y")
vlayer = QgsVectorLayer(uri, "asd", "delimitedtext")
QgsProject.instance().addMapLayer(vlayer)

and I got this in the console:

exec(open('/tmp/tmpu7edmqsu.py'.encode('utf-8')).read())

How to solve it?

PS:
I used QGIS 3.0.2 in Ubuntu 18. I also tried with the same QGIS version in Windows 8.1 and had the same results.


texts: Placing TIF image on top of WMS layer in QGIS

texts:
How to place TIF file on the top of WMS layer in QGIS?

When I import TIF file, it has different coordinate system, and it is placed far from WMS layer, but I want to put this TIF layer on the top of WMS layer.

I know how to make them transparent.


texts: Why is "RasterToOtherFormat" tool creating corrupted .TFW files?

texts:
Using the "RasterToOtherFormat" tool in ArcGIS 10.1 to convert tiff's with .tfwx file to tiff's with .tfw file, I get corrupt .twf files. That means, the saved coordinates in the new twf file are not correct and can't be read by QGIS.

Has this something to do with my ArcGIS Installation? If I run this tool from a different PC with ArcGIS 10.3 I get correct .tfw files. Is this a settings problem?


texts:
In QGIS 2.18 I'm try to sort features by an integer attribute with a Python script. The problem is, that some of the values of the sorting attribue are <= 0 (for sure they have to be <=0) and these features are ignorded in the sorting. How do I have to change the syntax to be sure, that also features with attribute [1] <= 0 are placed on the very first places of the table?

for feature in layer.getFeatures():
    count = layer.featureCount()
    geom = feature.geometry()
    attr = feature.attributes()

features = sorted(layer.getFeatures(),key = lambda feature: feature[1])
layer_provider = layer.dataProvider()
layer_provider.addFeatures(features)
layer.commitChanges()


texts: Creating a 3D voxel plot with LAD values in lidR

texts:
I am interested in the vertical distribution of the leaf area density (LAD) in my plot. That is why I am seeking to create voxels that represent LAD. I've been trying to achieve this with the grid_metrics3d() function. According to the documentation of this function the inputs necessary are a LAS object, a function to be applied to the voxels, and a resolution of the voxels.
The functions that can be applied to the voxels work the same way as the grid_metrics() function. One of the pre-made functions is the LAD function which computes a LAD value per voxel.

My limited experience in scripting got me this far:

library(lidR)
library(lidRplugins)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
voxelsLAD = grid_metrics3d(las, ~LAD(z=las@data$Z, dz = 1, k= 0.5, z0= 1), 5)
plot(voxelsLAD, color = "lad", colorPalette = heat.colors(5), bg = "black", legend = TRUE)

The resulting LAD values range from 0.00009805 to 0.247 which are in line with the values retrieved by the original research of which the function is based on. So far so good you would say, however plotting these values turned out to be difficult. The colorpalette does not apply itself to the voxels as can be seen in this image (all voxels retrieve the same color):

Right now I'm not sure whether this problem can be fixed by computing the LAD voxels otherwise (I might be misusing grid_metrics3d()), or by plotting the voxels otherwise.
FYI; I'm using Windows 10 with R 3.5.0


texts:
I am creating many models doing many tasks. Basically, the first model will join 2 tables and the steps are as follow: make table view for each table, then add join for, then it will copy the rows for the resulted table, then add a field, and then do field calculations and finally copy the resulted rows to produce the desired table.

Now one of the input tables has a fixed name, while the other could be named anything but it will always has the same fields, so only the name will change for that table. In the calculate field step, the issue I am having, is when I run the model using a new name for the input table; for example:
The first input table is named "ImpotTable" and this table always has the same name. The other input table (Join table) is named "DBO_Import" and it could be named "LPE_Import" or any other name. I want to do calculate field within the resulted join table for the field named "Number" which is from the join table.

So when I do calculate field, the "Field Calculator" window shows different field names under "Fields", so if I joined "DBO_Import" to "ImportTable" I will have "ImportTable.Number" field and "DBO_Import.Number" field; however, when I use the "LPE_Import" as a join table I will have "ImportTable.Number" and "LPE_Import.Number". So how I will tell the model to do calculate field for the "LPE_Import.Number" even if it was named "DBO_Import.Number" or any other join table name ("???.Number")?

After that, I am doing another model that is using the resulted join table as an input table. I am doing 4 different table selection, then copying the result from each step to finally produce 4 different tables based on the 4 different table select tools. Now the Issue I am having is similar to the previous one, I need to tell the model to perform the "Table select" on some fields even if the fields having different names.

For example one of the table select tools is based on the "Date" field, so if I used "DBO_Import" as a join field to produce the join table from the previous model, the "Date" field will be "DBO_Import.Date", but I want it to perform even if I used "LPE_Import" as a join table and the field is named "LPE_Import.Date", or having any name ("???.Date").

I am using ArcGIS 10.2.2 for Desktop.


texts: Device for Logging GPS and Accelerometer Data?

texts:
I am looking for a device that can track my location (lat,long) and speed all day long. Perhaps, I need GPS and accelerometer data bundled together.

What I want to do is to log my daily activities. I spend time walking, running, driving, taking a bus, staying in a building, etc., just like any normal human-being :) I am willing to carry the device with me at all times and would like to see every place I go into and my speed at certain positions.

My requirements are that the location and speed data should be recorded at least every ~20 seconds. I would also like the device to timestamp (date+time) each record. It would be enough if timestamp is precise to nearest minute and the location can go off about a maximum of 10-15m. (30-50ft.)

What I want is a portable device so I can carry all day long for about 2-3 months. I will only use it during casual, daily activities.

I don't want real-time tracking. I can collect the data from a storage later on.


texts: How to access shapefile table in QGIS Modeler R scripts?

texts:
I need to develop .rsx scripts for QGIS 2.0.1 Modeler to do statistical elaborations on data tables of loaded vector layers.

In this case I am not interested in accessing shapefile geometries, I'd like R to read layer table as a dataframe.

Unfortunately the only way (I know) to read data from the table of a vector layer in .rsx scripts is to state:

##vector_layer=vector

And then use the following R code

vector_df=as.data.frame(vector_later)

The OGR library used to access vector layers in R is very slow for shapefiles with many features.

Is there any other method to make R read the data table of a loaded vector layer in QGIS without accessing geometries?

Something like read.dbf() {foreign library} ....


texts: Get points from polygons with extreme angles

texts:
I intend to remove points from a polygons with extreme angles like the image link

Get the points
Remove from linestring

I make my buffer like this...

SELECT (ST_Dump(
          ST_Union(
             ST_BUFFER(the_geom, 10, 'endcap=flat join=round')))).geom as the_geom
FROM mytable

The error is on buffer unions, see the image.

I want to remove that points.


texts: Load an ArcGIS Vector Tiled Layer from a URL

texts:
.net core 3.1 Wpf application use Esri.ArcGISRuntime.WPF v100.7.0., tries to view vector tile from own tile server, but nothing displays.
my code:
ArcGISVectorTiledLayer vectorTiledLayer = new ArcGISVectorTiledLayer(uri);
var basemap = new Basemap(vectorTiledLayer);
Map map = new Map(basemap);
MapView.Map = map;

Similar code on javascript works fine.
Can you help me?


texts:
I'm currently doing a carte with CartoDB where I use CSS in order to show label but my issue is that I have entities which are very close and their label are difficult to read because of their superimposition.

I try some of the CSS properties' of CartoCSS (like text-min-padding or text-placement) but I failed.

Here is my CSS code

/** category visualization */

#doc {
   marker-opacity: 0.9;
   marker-line-color: #FFF;
   marker-line-width: 1;
   marker-line-opacity: 1;
   marker-placement: point;
   marker-type: ellipse;
   marker-width: 40;
   marker-allow-overlap: true;
}

#doc[type="ONG"] {
   marker-fill: #A6CEE3;
}
#doc[type="Think Tanks"] {
   marker-fill: #1F78B4;
}
#doc[type="Politique"] {
   marker-fill: #B2DF8A;
}

#doc::labels {
  text-name: [name];
  text-placement: point;
  text-face-name: 'DejaVu Sans Book';
  text-size: 20;
  text-fill: #000;
  text-halo-fill: #FFF;
  text-halo-radius: 1;
  text-dy: 0;
  text-allow-overlap: true;
  text-placement-type: simple;
  text-label-position-tolerance: 0;
}

Thanks a lot.
Have a nice day

A.H


texts: Clipping a shapefile in QGIS that has user defined projections

texts:
I have two shapefiles for Rwanda which have user defined projections. I am trying to clip the geology shapefile so that it has the same spatial extent as the districts shapefile. When I clip, I do get a saved file, but there are no data in it, the attribute table is empty as well.

Districts file has following:

  Generated CRS (+proj=tmerc +lat_0=0 +lon_0=30 +k=0.9999
  +x_0=500000 +y_0=5000000 +ellps=GRS80 +units=m +no_defs)

While geology file has the coordinate system:

  G * Generated CRS (+proj=tmerc +lat_0=0 +lon_0=30 +k=0.9996 +x_0=500000 +y_0=10000000 +ellps=clrk80 +units=m +no_defs)

I am having trouble to reproject them so that they gave the same coordinate system but if I do so, I am not able to view them together.


texts: I tried to import a dem.tif in a postgres database. I use postgres 9.5. First step was to create file via raster2pgsql. This file dem.sql (800.000kb), I'd like to import via psql into the database. The process start running, create sequence, and primary key, then CREATE TABLE...stop with ERROR out of memory character string with 859735225 Bytes to long for transformation ERROR Transaktion stopped etc. What can I do to bypass this problem?


texts: need to split some linestrings by points

texts:
please help me, i'm new in gis

1) i have а table with geom railways, which are multilinestrings (for example there are 2 ways one of which is interrupted, so this way will be multilinestrings consisting of many linestrings, wich have the same ID)

2) i have a table with points, which is on the line, but not in the linestring itself, i.e. each point is not a vertex of the linestrings

i need to split lines by points in this case (when 1 way is a multilinestring consisting of many linestrings, look at the picture below, column 'glput' is ID of the way)

an illustration of what I have and what i need

blue lines are segments that i must have

i try to use ST_LineSubstring (roads.geom, ST_LineLocatePoint (poits.geom1), ST_LineLocatePoint (points.geom2))

full query:

select ST_LineSubstring(st_geometryn, ST_LineLocatePoint (st_geometryn,geom1),
ST_LineLocatePoint (st_geometryn,geom2))
from table_road t1
join table_point t2
on t1.glput=t2.glput

but it is not working as it should, because ST_LineLocatePoint representing the location of the closest point on EACH of linestrings with the same ID, but it must interpolate point only to the line on which they lie


texts: How do I change the PostgreSQL data directory in the OpenGeo Suite 4.5?

texts:
I am using OpenGeo Suite and I would like to move the PostgreSQL databases files to a different drive than my operating system.  The default location installed by OpenGeo Suite is in the "C:\ProgramData\" directory, but I would like for it to be "F:\db\".

Some websites suggest that the data directory location can be changed by adding -D "F:\db\" to the ImagePath value of the PostgreSQL service key in the registry. However, for me this did not change where PostgreSQL is looking, and I have restarted the service several times.  Since OpenGeo customized the service and registry key names, I am afraid that they may have also customized the PostgreSQL service itself to prevent the changing of the data directory.

Does anyone have any knowledge of this?


texts: How do I load a UI file from a postgres database in QGIS?

texts:
I have a custom form I made and associated python code I made. Right now both are stored on mapped network drive and the layer style that gets saved to postgres references the .ui and .py file on the network drive.

Ideally, I would like it all to load as part of the style saved in the database. I was able to save the python code as part of the style by saving it in the "Provide code in this dialog" option but I cannot get a UI to load from the database.

When saving the style to the database, I did check the save UI to database and pointed it to the right file to load. But when I load the style, it still points to the same place on the network drive in the "Edit UI" input box.

I have tried various combinations like leaving it to auto-generate and whatnot but it doesn't seem to ever pull the form from the database no matter what I do.

Anyone know how to get this to work?


texts: Merge several polygons/features into one continuous area not working

texts:
I am using QGIS version 3.10. I have a farm boundaries layer (https://drive.google.com/drive/folders/1xYZmSboM75GACCWQyx-06wZ3EdHBcMrM?usp=sharing), and want to merge certain adjacent polygons into one fenceless area (the resultant merged area will be a nature reserve in which the fences between farms have been removed).

I tried to use the following method to achieve this: toggle editing -> clicked on polygons I want to merge -> merge selected features. A table "Merge Feature Attributes" pops up and from here I am not sure what to do, because, whichever option in the merge features table I select, the result is not what I want. For example, if I say skip all fields the result is that only some line segments have been removed:

The ultimate aim is to have the farm boundaries layer with the line segments between the selected polygons removed. I am analysing the movement of wildlife (which in the nature reserve is unobstructed by fences).
Additionally, as a second question, I would like to find out how to use a different method to achieve the same goal:
I have the same farm boundaries layer (green) and a national park layer (pink). I deleted the farms that belong to the nature reserve (white area):

Is it possible to save out the white area as a new polygon (without having to use digitising tools) and then merge it with the farm boundaries layer?


texts:
I have a basic ASP.net web app that loads an ESRI base map. The index page looks like this:

The JavaScript code:

      function (Map, FeatureLayer, LabelLayer, InfoTemplate) {
          var map;
          var initMap = function () {
              map = new Map("mapDiv", {
                  basemap: "topo",
                  center: [-82.6389, 27.872461], // longitude, latitude
                  zoom: 10
              });

              map.on("load", initMapLayers);
          };

          var initMapLayers = function () {

              var popup = new InfoTemplate();
              popup.setTitle("
NTA Median Household Income
");
              popup.setContent("Demographic Value Will Go here");
              var ntaFeatureLayer = new FeatureLayer("some/MapServer/0", {
                  mode: FeatureLayer.MODE_SNAPSHOT,
                  infoTemplate: popup,
                  outFields: ["NTA", "FranNTA"],
                  visible: false,
                  id: "coveredNTAs"
              });

              var ntaLabelLayer = new LabelLayer({ id: "ntaLabels", mode: "dynamic" });
              var featureLabelRenderer;
              ntaLabelLayer.addFeatureLayer(ntaFeatureLayer, featureLabelRenderer, "{NTA}");
              ntaLabelLayer.styling = true;

              var whereClause = "FranNTA IN ('3261LI','3261SV','3261SW','3261LJ','3261LD','3261LF','3261SU')";
              ntaFeatureLayer.setDefinitionExpression(whereClause);
              map.addLayer(ntaFeatureLayer, 1);
              map.addLayer(ntaLabelLayer, 4);
              ntaFeatureLayer.show();
              ntaLabelLayer.show();
              ntaFeatureLayer.setOpacity(0.5);
          };

          return {
              initMap: initMap
          }
      }
  );

The map loads up just fine but its height is prefixed to 400px.

Do i need a Custom CSS to fix this? Changing the styling in developer tools doesn't seem to change it.


texts: Generate annotation for features in ArcObjects

texts: I need to programmatically generate annotations for features.  I don't need the actual annotation feature, just the text placement, size, angle, etc. properties that are generated.  Those attributes are going to be stored in a different file format.  Is there a way to generate anno on the fly in ArcObjects to get at these properties?


texts: Saving and exporting filters in QGIS to reuse, is that possible?

texts:
I have a project with about 40 layers imported from .csv files. Each layer is data from Open Data Kit (ODK) surveys in a specific geographic area. There are about 15 filters that people want to apply to each layer. They don't need to apply them to all layers at once. What they want to do is look at a layer, apply a filter, then perhaps compare to another layer. The filters themselves won't change. It can actually become quite tedious to add a filter, copy and paste the filter and then do same for a whole bunch of layers.

To me, an ideal situation would be to apply a filter, export the filter as a text file perhaps, and then use an import that filter, similar to what you can do with a specific style. I don't know if this is possible in QGIS.

I took a look at this question and it helps, sort of gets to what I want to do, but rather than apply a filter to multiple layers, I would like to be able to select a pre-coded filter.


texts: Add each Landsat image as a seperate image band in Google Earth Engine

texts:
I would like to stack all Landsat 8 images from a certain time frame in Google Earth Engine by adding the blue band from each image as a separate band and exporting as a single image.

How would I go about this? I am able to load the collection in, but am unsure what to do next.

// Load Landsat 8 ImageCollection for area.
var collection=ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .filterDate('2015-01-01','2015-12-31')
  .filterBounds(AOI).select('B2');
print('Landsat: ', collection);


texts: cpu usage goes up when i move mouse to show xy coordinates

texts:
I want to show x , y coordinates when i mouse move over map in a webgis .The problem is whenever i move the mouse, the cpu usage goes up to a maximum of 40%.Is there any solution to reduce cpu usage?

I'm using ArcGIS API silverlight 3.1 , ArcGIS server 10.1 , sql server 2008 R2 , windows 7


texts: Splitting shapefile into separate files within Meridians using QGIS

texts:
Given a .shp that covers the globe, how could QGIS be used to segment that into 90 .shp files, each file containing lines/points/polygons/multipolygons that are between the first and second then second and third then third and fourth, etc?
Issue is a series of gigantic source files with global coverage that have no common attributes that would allow them to be conventionally partitioned. Since the source files have a relatively even distribution (not over the Pacific ocean of course), cutting them up into smaller files will make them cleaner to process.
I was hoping there was a hidden command in ogr2ogr that would do this, but I don't see one.


texts: Tilemill: Render a raster only where a shape layer is underneath

texts:
I'd like to mask a raster to only those areas, where there is a shape layer underneath. To be more precise: I have a land layer and a hillshading on top and I only want my hillshading to appear on top of the land, not it the water areas.

I tried around with different comp-ops but i didn't get what I want.

The result I'd like is see the land shape (it's color) underneath and the raster only visible on top where there is a land layer underneath.

Is that somehow possible?


texts: Combining different geometry types into same feature class?

texts:
I've searched the question in this blog. Maybe I don't use the right words but anyway my problem is not present in any post so I'll try to explain my question by opening a new one.

I have a polygon shape (eg states).
I calculated the centroids and the function saved them in an another file (Feature To Point function).

Now I have 2 separated layers, point and polygon shape (centroids and polygons), but I need a unique shape.
Functions Merge and Union (if I understand what I've read) can only be made between shape with the same feature class (point-point, line-line, polygon-polygon).

I need a layer that contains, for example, the polygons of the states and their centroids
Is there a way to get the results I need?
Or a procedure to get the centroids of the polygons on the same shape?
I hope I explained
Sorry for the bad english!


texts:
I have two SHP files in QGIS project (points: NODES.shp and lines: BRANCH.shp). NODES.shp has one ID attribute NODE_CON. BRANCH.shp has attributes NODE1 and NODE2 which represent start and end node of topology network.
I need the tool wich at inserting feature into BRANCH.shp automaticly detect NODE_CON of the start node of the new polyline and write it into NODE1. Similar it need to detect NODE_CON of the end node and write it into NODE2.
My idea:
I am trying with Default value on BRANCH form. My idea is to write the function getNodeCon with three parameters:

layer_name,
field_id and
geometry of point.

I would call the function like that:
getNodeCon('NODES', 'NODE_CON', start_point($geometry)).
---------function definition start---------------
import os
from qgis.core import *
from qgis.gui import *

@qgsfunction(args=3, group='Custom')
def getNodeCon(layer_name, conn_field, tocka):

    layerList = QgsProject.instance().mapLayersByName(layer_name)
    layer = layerList[0]

    NodeCon = 'NOK'

    for feat in layer.getFeatures():
            if feat.geometry().within.tocka.geometry().buffer(0.5,100):
                    NodeCon = feat.fieldNameIndex(conn_field)
                    break
    return NodeCon

---------function definition end---------------
I get error measage:
Eval Error: QgsProject.mapLayersByName(): argument 1 has unexpected type 'list'


texts: Why krige function return all NA's values for the prediction when trying to interpolate a variable?

texts:
I am trying to build a map with kriging interpolation for a variable. When running the krige function the return is all NA values and the kriged map plot only one value. I couldn't find out why and I am a really basic user of these tools. Some older related answers mentioned the possibility that duplicate observations were causing this, but I checked and there are no duplicate records. I was able to build the variogram.
Here is the base map (shp file in a .zip file): http://www.filedropper.com/lpr000b16ae
Here is the data I want to interpolate: http://www.filedropper.com/yielddata
library(tidyverse)
library(sf)
library(gstat)
library(stars)

## a) Defining the projection of the maps
proj_lambert <- "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

## b) Map Canada selected provinces.
canada.map <- st_read("C:[your path to]/lpr_000b16a_e/lpr_000b16a_e.shp") %>%
  filter(PRNAME == "Manitoba" |
           PRNAME == "Alberta" |
           PRNAME == "Saskatchewan") %>%
  st_transform(proj_lambert )

## c) Data set
data.yield <- read_csv("C:[your path to]/yield.data.csv")
data.yield.coords <- data.yield %>% select(lat, long)
data.yield.data <- data.yield %>% select(av.yield)

#### Creating LAT LONG SpatialPoints
coordinates(data.yield.coords) = c("long", "lat")
proj4string(data.yield.coords) <- CRS("+init=epsg:4326") # 4326 is the EPSG identifier of WGS84.

### Converting the LAT LONG to the lambert CRS shown above
lambert_data.yield.coords = spTransform(data.yield.coords, proj_lambert)
lambert_data.yield.coords = as.data.frame(lambert_data.yield.coords)

#### Adding the Type column back to the data frame, with the new polar coordinates
lambert_data = cbind(lambert_data.yield.coords, data.yield.data)

data.yield.sf <- st_as_sf(data.yield, coords = c("long", "lat"), crs = proj_lambert)

## d) Variogram fit
v = variogram(av.yield~1, data.yield.sf)
plot(v, plot.numbers = TRUE)

v.m = fit.variogram(v, vgm("Sph", psill = 1, range = 5, nugget = 1))
plot(v, v.m, plot.numbers = FALSE)

This is the variogram:

And the rest of the code
## e) Building a empty grid
st_bbox(canada.map) %>%
  st_as_stars(dx = 1000) %>%
  st_set_crs(proj_lambert) %>%
  st_crop(canada.map) -> grd

## f) Kriging interpolation
k = krige(av.yield~1, data.yield.sf, grd, v.m)

## g) plotting all together
ggplot() +
  geom_stars(data = k, aes(fill = var1.pred, x = x, y = y)) +
  geom_sf(data = canada.map, color = "gray", fill="transparent") +
  geom_point(data = lambert_data, aes(x = long, y = lat), size = 1) +
  theme_bw()

And here is the final map


texts:
Why my mobile GPS routing calculates routes almost instantly and using pgRouting in my robust PC is so time costly ?

EDIT:

My test:

In my cellphone (Android 4.2.2 - 1G RAM - Dual core - Arm v9) I have iGo Primo.
In my PC ( 8G RAM - Core i5 ) I have Postgre + PostGIS + pgRouting

Table OSM Planet "planet_osm_line" with previous work:

alter table planet_osm_line add column source integer;
alter table planet_osm_line add column target integer;
select pgr_createTopology('planet_osm_line', 0.0001, 'way', 'osm_id');

select * from planet_osm_line where name like '%BR-040%' // (source : 2751784)
select * from planet_osm_line where name like '%BR-101%' // (target : 2040999)

SELECT seq, id1 AS node, id2 AS edge, cost, name, way
  FROM pgr_dijkstra('SELECT osm_id as id, source, target, st_length(way) as cost FROM planet_osm_line',  2751784, 2040999, false, false  ) as di
  JOIN planet_osm_line pt
  ON di.id2 = pt.osm_id ;

The route is (both Brazil Highway):

Source: BR-040

Target: BR-101

My cellphone do the job in less than two seconds and my PC takes more than 3 minutes (only the SELECT route - pgr_dijkstra - part).


texts: PostGIS - spatial index not used on ST_Intersects when doing update

texts:
I'm trying to update the sites table with any ancient woodlands that each site intersects. I'm running the SQL that looks like below, which works but is very slow and doesn't appear to use spatial indexes. Is there a way to rewrite it that would use them?

UPDATE sites
SET woodlands = (
    SELECT name
    FROM ancientwoodlands, sites
    WHERE ST_Intersects(
        (SELECT sites.geom WHERE sites.key_column = 1),
        ancientwoodlands.geom
    )
    LIMIT 1
)
WHERE key_column = 1;


texts: What analysis tool I can use to identify the relation between points and areas

texts:
I want a binary result for each purple point in this map, i.e. whether in or out od the blue polygon areas?

In this case, one is in and the other is out. What analysis tool in ArcGIS Pro (also arcpy package) I can use for this purpose?

I try to use AggregatePoints_sfa as follows:

AggregatePoints_sfa('Point_Location','Areas','Result')

But an error was incurred:

OSError: The toolbox file https://analysis2.arcgis.com/arcgis;services/tasks;UseSSOIdentityIfPortalOwned was not found.

Not sure if it is something to do with my ArcGis PRO version which is 2.2.0.


texts: Downloading Sentinel-2 data from shapefile of whole country

texts:
I have shapefile of a small country: Portugal
I need to download Sentinel-2 data corresponding to the shapefile, or download a bigger one and then cut it.
I tried using SCP plugin of QGIS but only know how to download a small rectangle


texts: Coordinates system in degrees in plugins for QGIS

texts: After adding a vector layer of coordinates of settlements of a part of the world INTO QGIS 1.8, I wish to add a topographic or other nice looking, colored map. Trying plugins like maps from Google Earth changes the coordinate system into meters, This has the consequence that I cannot lable coordinates in degrees with a grid. Are there plugins of a world map that can be displayed in degrees in QGIS?


texts:
I currently have a shp file with ~300 points. This shp file stores a bore name, x, y, z, and a few other attributes.

In a test PostGIS database, via QGIS, I have imported that shp file, called 'collars'. I have then gone and added holes to the 'collars' table, adding 3 new rows with all the required data, boreid, x, y, z etc.

Unfortunately I have realised upon re-opening QGIS that the layer does not automatically display these new holes. I understand that as I didn't add anything into the 'geometry' column these news rows of data have no spatial locations in the PostGIS database.

What I would like to know is can I cause the layers spatial data to recalculate, based off the x and y values, every time I open QGIS or whenever new data is added so that this 'collars' layer becomes a 'live' linked layer displaying new data as it is added in?


texts: Using selection of features in ArcMap in Python script?

texts:
I am trying to write a python script and am stuck on the first obstacle.

I want the user to select features in ArcMap and then use these selected features in a python script.

Is that at all possible and does someone have a code example?

The code that I have now uses all the features instead of just the selected ones. Please see code snippet below.

arcpy.MakeFeatureLayer_management(pipes,"pipeslyr")
arcpy.SetParameter(0,"pipeslyr")
pipesnew = os.path.join(ws,"pipesnew")

arcpy.CopyFeatures_management("pipeslyr", pipesnew)


texts: Extract pixel value for corresponding gps location in R

texts:
I have animal gps data with several locations for each day and regular gaps of several days when no location of the animal was recorded. Additionally, I have satellite data in intervals of 16 days. I would now like to extract the pixel value that corresponds to the specific point and to the specific time.

This means if the location of the animals was recorded e.g. 2 days before a satellite image was taken, I want to extract the pixel value of this image (that comes after) rather than from the images that was taken 14 days before the animal location was recorded. I always want to extract from the image that lies closer according to the time.

I created some test data that hopefully illustrates the problem:

library(sp)
library(raster)

### Create test data
# create first raster
edc2012001_m <- raster(ncol=36, nrow=18)
edc2012001_m[] <- sample(1:ncell(edc2012001_m))

# create second raster
edc2012017_m <- raster(ncol=36, nrow=18)
edc2012017_m[] <- sample(1:ncell(edc2012017_m))

rasters<-stack(edc2012001_m,edc2012017_m)

# Create xy coordinates
time<-c("2012-01-01", "2012-01-01", "2012-01-01", "2012-01-02", "2012-01-02", "2012-01-02", "2012-01-12", "2012-01-12", "2012-01-13", "2012-01-13")
x <- rep(-50,10)
y <- sample(c(-80:80), 10)
data<-data.frame(x,y,time)

# Convert data to spatial points data frame
coordinates(data) <- c("x","y")

### Extract all data from raster stack
extract(rasters, data)

The number string in the name of the first raster indicates that the image was taken on the first day in 2012, the second image was taken on the 17th day of this year.

The 7th location in the test data for example shall now be extracted from the second raster file as it lies closer according to the time.

All in all I have 87 raster files and 600 observations. I am really not sure how to program this. I guess I could use substr() to retrieve the date info from the raster names. But besides that...
I am thankful for every hint I can get, also about functions that might be helpful in this context.


texts: Passing extent from web map to PDF for integration?

texts:
I'm trying to figure out whether I can take an extent (or coordinatess with a Zoom level included) to a PDF in order to author a PDF that displays both GIS and report-styled business data?

Is this something possible with ArcGIS API for JavaScript or even via Flex?


texts:
I have some doubts with the generation of a cost surface with r.walk of GRASS. It's assumed that the resulting layer is measured in seconds, but when I try to convert seconds in days march (6h per day), the resulting layer has no sense.

This is the workflow I'm following:

1st.- I have several vector layers as roads, rivers, bridges, and a raster layer of elevation.

2nd.- Rasterize the vector layers with a field named "weight" (PESO in spanish) where I have weighted this variables according to if they facilitate the walking or not. Reclassify is made if it's needed.

3rd.- Find Slope (in percent) layer and reclasify as follow:

0-10=1;
11-maxValue=100

4th.- r.cross with all the rasterized layer. (roads, bridge, river and reclassified slope).

5th.- I create a new friction layer using this formula proposed by Pandolf (1977) with r.mapcalc :

M=1.5W+2.0(W+L)(L/W)^2+N(W+L)(1.5V^2+0.35V*abs(G+10))

6th.- With this layer measured in Watts I use r.walk with default parameters.

7th.- Once obtained this layer measured in seconds I try to convert it in hours dividing 3600 and then convert these hours in days march dividing 6 (6 hours per day).

Well, the results are not what I expected, the max data of this layer is about 8.000... No sense, it should be something like 3 or 4 days, maybe.

Here are the files I'm using to do this workflow.

I don't know what I'm missing, maybe is the 4th step wrong? Maybe the 5th?

Another point of view would be nice to me. If you need more info, please ask for it and I will add to the question.

Thanks a lot!


texts: ogr2ogr get geojson from map server misses features and creates duplicates

texts:
I'm querying this map service using ogr2ogr to create a local copy of it's features:
http://geoportal.menlhk.go.id/arcgis/rest/services/KLHK/Deforestasi_2006_2009/MapServer/0

Here's my ogr2ogr code, using where clause objectid > 0 to select everything:

ogr2ogr 2006_2009.geojson -f GEOJson "http://geoportal.menlhk.go.id/arcgis/rest/services/KLHK/Deforestasi_2006_2009/MapServer/0/query?where=objectid>0&outfields=*&f=json" OGRGeoJSON -gt 1000 -s_srs EPSG:3857 -t_srs EPSG:4326

My output (2006_2009.geojson) has 8724 features, which matches the record count here:
http://geoportal.menlhk.go.id/arcgis/rest/services/KLHK/Deforestasi_2006_2009/MapServer/0/query?where=1%3D1&returnCountOnly=true&f=json

However, when I compare the shapefile in Arc with the shapefile in ArcOnline (coming directly from the source map server) I found the shapefile I created was missing features.

The reason I was seeing matching feature counts is because the ogr2ogr output has duplicate records:

ogrinfo 2006_2009.geojson OGRGeoJSON -dialect sqlite -sql "SELECT objectid, count(objectid) FROM OGRGeoJSON GROUP BY objectid HAVING (COUNT(objectid) > 1)"

The above query identifies 147 instances of duplicate objectids, suggesting that we're missing 147 records of real data given that the record counts match.
Looking  at the output dataset, I can see that there's nothing for objectids between 852 and 892. I can use ogr2ogr to extract these features successfully:

ogr2ogr -f GeoJSON missing_features.geojson "http://geoportal.menlhk.go.id/arcgis/rest/services/KLHK/Deforestasi_2006_2009/MapServer/0/query?where=objectid%3E852%20AND%20objectid%3C892&outFields=*&f=json" OGRGeoJSON -gt 1000 -s_srs EPSG:3857 -t_srs EPSG:4326

But for some reason these are not present in the original output, where the where clause selected all features with objectids > 0.

In sum, two issues:

Why doesn't ogr2ogr export objectids 852 - 892 (and others) given a where clause of objectid > 0?
More troubling-- why does it duplicate features to make the record counts equal?


texts: I'm trying to create a join where I take the sum of the values of one of the attributes in the point layer within each buffer and add it to the bufferlayers attribute table. I can't figure out which settings I need. Do I get all values with intersect or should i choose contains or touches? I want the aggregated value from one of the attributes within so I suppose I should take the "Attribute summary" that you can se in the picture?