texts: Continuous raster legend in QGIS map composer

texts: Is it possible to make continuous raster legend in QGIS? Something like in ArcGIS???


texts: Distance in polygon from points in postgis?

texts: I've got two GeoJSON files. One is points (gas stations) and the other is the local administrative area. I'm wondering is it possible to calculate the farthest distance in that polygon from the point?


texts: Adding new menu item to QGIS Desktop app?

texts:
I need to add new menu item to the top level menu of QGIS Desktop (2.7.0);
I.e. it shall be on same level as "Project", "Edit", "View", "Layer" etc.

As far as I know there is pretty nice way how hide menu items (Settings->Customization)
Plus using QgisInterface (http://qgis.org/api/classQgisInterface.html) I can add new items to Menus already defined in QGIS (addLayerMenu, addPluginToDatabaseMenu etc) from my plugin code.

But I need new Menu item on top layer (this is req from customer).


texts: ENVI Example Based FX segment selection slow

texts:
I am using the Example Based FX workflow in ENVI 5.5 and the segment selection tool is extremely slow. It used to be fairly fast but now takes a long time (30 seconds or more) to process one small segment even one as small as a 10m pixel.
Does anyone know of a setting to speed this up or any tips for correcting this slow behavior?


texts: How to convert decimal point to comma in QGIS?

texts: I am working with data in Colombia where decimal points are written as commas. I would love to create a new column in QGIS that would convert all decimal points (periods) to commas for the purpose of labeling. Any thoughts?


texts: How to export a PDF from composer when I use the Google Physical layer as a basemap in QGIS

texts: I am unable to open a new map I've just made in print composer if I use the Google Physical basemap layer. After I export the map I've made and try to open it on my desktop, it reads that the PDF is damaged and cannot be opened. I tried using other topographic base maps from the open layers plugin as well and had similar problems. The export works perfectly fine with arial and hybrid layers. Any suggestions?


texts: Determining when layer visibility is turned ON or OFF (changed) in QGIS plugin?

texts:
I need to update a vector layer's data when the visibility is turned back on.  My update does not occur when the visibility is off.

To what signal do I need to connect to catch the visibility change?


texts: How to place points along a line in a specific offset using python / arcpy?

texts:
I have a street shapefile with information about the street name, house number from / to on the left side and house number from / to on the right side.

My aim is to create a point shapefile that visualizes all houses along the streets. A first prototype is already implemented. With a Search-cursor I am reading out the total number of houses by building the range by using the paramaters house number from / to on the left side and house number from / to on the right side. The total amount is divided by the length of the street so that I have a value that show me in which distance I have to place the points along the line to map all houses.

The result is that all house numbers are placed like a string of pearls.

Is there a way to place the numbers with an specific offset along the line?

The script is written in python using arcpy and ArcGIS 10.1 SP1.

Thanks!


texts: How to calculate distance on large raster in R?

texts:
I'm trying to do some raster algebra and one of my layers is a distance from a coastline layer. I'm using rasters with 12.5m resolution and when I use the distance function, I seem to run into memory issues. Is there a way to set a maximum distance like dist in gDistance? Are there other ways to speed this calculation up? Here's a mock raster with similar specs:

#build blank raster modeled after large raster
r <- raster(xmn=1035792, xmx= 1116792, ymn=825303.6, ymx=937803.6, resolution = 12.5,crs = "+init=epsg:3174")
r <- setValues(r, 0)

#create lake polygon
x <- c(1199999, 1080000, 1093067, 1090190, 1087977, 1070419, 1180419)
y <- c(957803.6,937803.6, 894366.9, 872153.9, 853703.0, 825353.6, 805353.6)
poly.lake <- SpatialPolygons(list(Polygons(list(Polygon(data.frame(x,y))), ID = 1)))

# make values NA where lake polygon does not intersect raster
r <- mask(r, poly.lake)

#run distance function
r.dist <-distance(r)


texts: Bulk field renaming in file geodatabase?

texts:
I have a file geodatabase with several hundred feature classes in it.

Some feature classes have a field that got named incorrectly, and it got named "SCHEMA_FEATURECLASS_ENTITY" instead of just ENTITY.  Consequently, their final destination (Oracle) complains about the very long field names.

I'd like to rename those fields to have their proper short name.  I did find Changing feature class and field aliases in bulk using ArcPy? but (from what I understand) renaming a field is an add-copy-delete operation rather than just setting a property.

Has anyone got a quick-and-dirty method for doing a bulk rename?  I have ArcCatalog 9.3.1 but NOT Visual Studio (client environment...).


texts: Cleaning large Shapefile using v.clean in order to dissolve features?

texts:
I'm having real problems cleaning up and reducing the size of a Shapefile which is published by the UK's Environment Agency.

It shows the extents of the LiDAR open data they publish: each polygon is a survey flight, and there are fields for the resolution, the date of the flight and so on.

Eventually I'd like to end up with a dissolved (merging all features) and simplified Shapefile which is much smaller in size. At the moment:

The .shp is 360.6MB
The .dbp is 26.9MB
There are 121,753 polygons

I think one of the reasons the file is so large is that there are many small 'specks' of data (which aren't important for my purposes):

What I've tried so far:

Dissolving with QGIS: this seemed to make no progress on such a big file so I cancelled it after a while
Dissolving with OGR (ogr2ogr dissolved.shp LIDAR_composite_extents_2015.shp -dialect sqlite -sql "SELECT ST_Union(geometry) AS geometry FROM LIDAR_composite_extents_2015"): after a few hours I would get an error like this: GEOS error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 221912.50000000093 50580.000000001863 at 221912.50000000093 50580.000000001863 – I guess this will just be the first of many errors
Cleaning the Shapefile, first with QGIS's Check Geometry Validity, then with GRASS's v.clean (I tried bpol), but the cleaned file still fails dissolve (I also tried adding a zero buffer)
Converting the multipart polygons to singleparts, adding a geometry column and removing features smaller than a certain area. This took about 100MB off the filesize, but didn't affect the 'holes' within polygons – to get rid of those I tried to make a difference layer, but the difference operation consistently fails (if I check ignore invalid geometry, it creates lots of features which are visible in the attribute table, but doesn't display them).

I will simplify the file in the end, but don't want to do this before dissolving in case I introduce slivers.

Should I use a different tool within v.clean?

I have not yet tried to split the file into regions and dissolve those, then dissolve the regions.


texts: How to sample points on specific landcover class using Image.stratifiedSample in Google Earth Engine?

texts:
I need to sample NDVI time series on the crop classe defined in the Copernicus GLDAS Landover dataset (value: 40). The logical approach seems to be using Image.stratifiedSample on the landcover, retrieve a featureCollection and use it to sample NDVI from another dataset. However I get a EEException: An internal error has occurred upon defining the stratified sample with the code below. Can anyone spot what is wrong? Or is there a better approach I haven't seen?
import ee
ee.Initialize()

xmin = -72
xmax = -67
ymin = -42
ymax = -39

ROI = ee.Geometry.Rectangle([xmin, ymin, xmax, ymax])

coper = ee.ImageCollection("COPERNICUS/Landcover/100m/Proba-V-C3/Global").first().select('discrete_classification')
coper = coper.clip(ROI)

points = coper.stratifiedSample(
  numPoints=10,
  classBand='discrete_classification',
  region=ROI,
  scale=100,
  classValues=[40]
)

task = ee.batch.Export.table.toDrive(collection=points, fileFormat='GeoJSON')
task.start()


texts: Extract raster information from a spdf and for different buffer size in a loop

texts:
For a project in geoprocessing, I need to extract information from a raster but for different buffer sizes. I have a raster of the population (100 m size pixel) and a SpatialPolygonDataFrame with the localisation of volcanoes. According to the volcanic explosivity index (VEI) of those volcanoes, I created a new column in the spdf with different buffer sizes (VEI = 1 --> buffer_size = 5000; VEI = 2 --> buffer_size = 10000, ...). Now, I want to extract the sum of pixel included inside the buffer. I want to do it the more efficient ! Then, I wanted to make a loop, but it doesn't work... Can somebody tell me why?

When I want to extract for a single buffer value, there is no problem (see the following code) :

people_number <- extract(x = population_proj, y = locations_spdf_proj, buffer = 10000, fun = sum)

But When I want to do that in a loop, I have a warning message (Error in .cellValues(x, y, ...)) :

  for (i in 1:length(locations_spdf_proj)){

  people[[i]] <- extract(x = population_proj,                            # raster layer
                         y = locations_spdf_proj[[i]],                   # SPDF with locations for buffer
                         buffer = locations_spdf_proj$buffer_size[[i]],  # buffer size
                         fun = sum)                                       # value to extract (sum of all pixels)

}


texts:
I tried to do this code Array to Raster image but there was the following error

ERROR 4: Unable to open myraster.tif to obtain file list.
ERROR 4: Unable to open EPSG support file gcs.csv.
Try setting the GDAL_DATA environment variable to point to the
directory containing EPSG csv files.

this is the code

import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

array = np.array(( (0.1, 0.2, 0.3, 0.4),
                   (0.2, 0.3, 0.4, 0.5),
                   (0.3, 0.4, 0.5, 0.6),
                   (0.4, 0.5, 0.6, 0.7),
                   (0.5, 0.6, 0.7, 0.8) ))
# My image array
lat = np.array(( (10.0, 10.0, 10.0, 10.0),
                 ( 9.5,  9.5,  9.5,  9.5),
                 ( 9.0,  9.0,  9.0,  9.0),
                 ( 8.5,  8.5,  8.5,  8.5),
                 ( 8.0,  8.0,  8.0,  8.0) ))
lon = np.array(( (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5) ))
# For each pixel I know it's latitude and longitude.
# As you'll see below you only really need the coordinates of
# one corner, and the resolution of the file.

xmin,ymin,xmax,ymax = [lon.min(),lat.min(),lon.max(),lat.max()]
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
# That's (top left x, w-e pixel resolution, rotation (0 if North is up),
#         top left y, rotation (0 if North is up), n-s pixel resolution)
# I don't know why rotation is in twice???

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.
                                             # Anyone know how to specify the
                                             # IAU2000:49900 Mars encoding?
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system
                                                   # to the file
output_raster.GetRasterBand(1).WriteArray(array)   # Writes my array to the raster


texts: Rounding to nearest multiple of pixel resolution using ArcPy

texts:
ArcPy can return the pixel resolution:

arcpy.GetRasterProperties_management(rasterfile, "CELLSIZEX", "").getOutput(0))

Let's say this returns "3.2808"

I wish to be able to take a given input, say "10" for example, and return the nearest multiple of the pixel resolution

n = 3.2808
MagicRoundingFunction(n, 10)
>>> 9.8424

Is there a built in function to accomplish this or must a custom function be developed?

If one exists, I would rather save time using something built in.


texts: Passing selected polygon to INTERSECTS filter in a WFS request

texts:
In the code below,

   <ogc:Filter>
     <ogc:Intersects>
       <ogc:PropertyName>the_geom</ogc:PropertyName>
       <ogc:Function name="querySingle">
          <ogc:Literal>sf:restricted</ogc:Literal>
          <ogc:Literal>the_geom</ogc:Literal>
          <ogc:Literal>cat = 3</ogc:Literal>
       </ogc:Function>
     </ogc:Intersects>
   </ogc:Filter>

I'd like to replace

   <ogc:Literal>sf:restricted</ogc:Literal>
   <ogc:Literal>the_geom</ogc:Literal>
   <ogc:Literal>cat = 3</ogc:Literal>

with a polygon selected by a user. I get the polygon from a selection event. I've tried a number of options, but none is working, yet. Any assistance will be highly appreciated.

My environment: GeoServer 2.1.3, OpenLayers 2.11.

I've transformed my geometry to EPSG:4326, the one used by my layers in GeoServer:

var geometry = event.feature.geometry.transform(new OpenLayers.Projection("900913"), new OpenLayers.Projection("EPSG:4326"));

The event.feature.geometry represents a polygon selected by a user.
I've also converted the geometry to WKT format as shown below:

var wktParser = new OpenLayers.Format.WKT();
var feature = new OpenLayers.Feature.Vector();
feature.geometry = geometry;

var wktGeom = wktParser.write(feature);

I've then passed the geometry in WKT format to the filter as shown below:

'<ogc:Filter>' +
'<ogc:Intersects>' +
    '<ogc:PropertyName>the_geom</ogc:PropertyName>' +
    '<ogc:Function name="querySingle">' +
        '<ogc:Literal>ws:layer_2</ogc:Literal>' +
        '<ogc:Literal>the_geom</ogc:Literal>' +
        '<ogc:Literal>INTERSECTS(the_geom,' + wktGeom + ')</ogc:Literal>' +
    '</ogc:Function>' +
'</ogc:Intersects>' +
'</ogc:Filter>';

But when I send a request, the following error is returned:

<ows:ExceptionText>Process failed during execution
java.lang.ClassCastException: org.geotools.filter.AttributeExpressionImpl cannot be cast to org.opengis.filter.Filter
org.geotools.filter.AttributeExpressionImpl cannot be cast to org.opengis.filter.Filter</ows:ExceptionText>.

Any idea on what I may be doing wrong?


texts: Getting new start and end points (X & Y) for new clipped line sections using ArcGIS Desktop?

texts:
I'm trying to get start and end points, and lengths for clipped sections of lines (roads) that fall in my polygon. Each road has two to three different sections within the polygon. I clipped the roads with the polygon, to show me only sections of the roads that go through the polygon.

Now I am trying to break up each section into its own line segment.

Then I am trying to get the start and end points (X&Y coordinates), and length of each segmented line.

I tried linear referencing, but has not work, might be doing it wrong. Also each road has a start point and a end point, but those are for the entire road, these points are not within my clipped area.

I'm using ArcGIS Desktop 10.2.


texts:
I'm trying to build an interactive game map with Leaflet and I would like to implement nested layer groups. I tried doing it the obvious way: putting some layer groups inside of a "master" layer group, and that does kinda work, but the control won't work with it even if I write custom event handlers.
I spent the past hour looking into the source of Leaflet trying to understand why this won't work, but I can't figure it out.

// Load markers
const overlay_materials = {};
const materials_lg = L.layerGroup();
overlay_materials["Materials"] = materials_lg;

for(let key in markers['materials']) {
    let material = markers['materials'][key];
    let layerGroup = L.layerGroup();
    for(let coords in material['coords']) {
        coords = material['coords'][coords];
        layerGroup.addLayer(
            L.marker(coords, {
                icon: L.icon({
                    iconUrl: markers['icon_base'] + material['icon'],
                    iconSize: material['size']
                }),
                title: key,
                alt: key
            }));
    }

    materials_lg.addLayer(layerGroup);
    overlay_materials[key] = layerGroup;
}

materials_lg.addTo(map);

// Add control for layers
L.control.layers(null, overlay_materials, {
    collapsed: false
}).addTo(map);

materials_lg.on('add', e => {
    materials_lg.eachLayer(lg => {
        map.addLayer(lg);
    })
});

materials_lg.on('remove', e => {
    console.log("Removed");
    materials_lg.eachLayer(lg => {
        console.log(lg);
        map.removeLayer(lg);
    })
});

If all the sub layers groups are not selected it does work fine, but it does not update the control checkboxes, so if one of the other layer groups is selected it won't work as expected.


texts: List of the United States (US) state boundaries / borders as latitude/longitude pairs for geofence?

texts:
I'm just trying to find a single simple file (CSV or XML ideally) which has the latitude and longitude pairings for each of the 50 United States (plus DC, though the addition of provinces such as Puerto Rico and Guam, or even Canada and Mexico wouldn't be unwelcome).

Is there such a list available somewhere?  I've tried searching and browsing OpenStreetMaps, downloading the US.zip from GeoNames, the Census website, etc., but everything seems to come in shape files and such which I'm not terribly familiar with.

My end goal (if it helps) is to simply look-up the state by latitude/longitude using my LAMP stack (Linux Apache (actually using Nginx) MySQL and PHP), though I will have 1000s of records per hour, so using an external web service isn't feasible.


texts: Default PostGIS type of geometry in SELECT

texts:
Is there an option in PostGIS (in query or table option) to set the return type of a Geometry?

I'm using this query:

select * from "locations" where "locations"."id" = 2 limit 1

And that will return

2;"My Home";"My Address";"0101000020E6100000EC14AB0661855EC010CF126404B64240"

But I want to default return the geometry as JSON, without naming the parameters, I mean, I don't want to do:

select id, name, address, ST_AsGeoJSON (location) from "locations" where "locations"."id" = 2 limit 1

I just want to continue using select * and returning location as GeoJSON by default. Is this possible?


texts: Misalignment between overlay and base layers in Leaflet web map

texts:
I use esri-leaflet (https://esri.github.io/esri-leaflet/api-reference/layers/dynamic-map-layer.html) library for LEaflet to display ESRI Dynamic Map Layer (http://gismaps.baladia.gov.kw/arcgis/rest/services/KA/ParcelLocation_MKM/MapServer) which have not EPSG:3857 coordinate system (PROJCS["KTM_KM",GEOGCS["GCS_KUDAMS",DATUM["D_Kuwait_Utility",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",499975.0],PARAMETER["False_Northing",209269.0],PARAMETER["Central_Meridian",48.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Central_Parallel",29.0],UNIT["Meter",1.0]])

It is possible to display this layer on Leaflet web map without misalignment?

Web map example https://jsfiddle.net/8y07mcvb/


texts:
I have a Dockable window control that I want to have several properties associated with it. Specifically a data source. I currently have to hard code an Image service or a file GDB at a specific path. I want to create a option that allows a user to choose between an Image Service via URL or a file GDB given a local .gdb file.

I found Property Page, and figured it would be my best best, but I feel like I'm missing something in the documentation that describes how to actually use it.

I've reviewed Creating property pages, and How to create property pages and property sheets and I've been able to open a property page with my custom property sheet, but now I dont understand how I'm supposed to 'use' the property page.

How do I properly utilize the PropertyPage to update custom properties in a custom Desktop Add in, or what else should I be using to accomplish this in ArcObjects?

I've considered having the options to change the data source on the Dockable window itself, but it already has a lot of controls and dont want to add confusion.


texts: Increase contrast of merged Planet Labs scenes

texts:
I've download a few PlanetScope scenes from Celaya, MX.

Taken alone they have a good contrast.

Then I've merged them with gdalbuildvrt, which lowered the contrast: Link

I've thought, that I could revert this again by stretching the contrast after this guide. gdalinfo -mm test.vrt gave me:

Band 1 Block=128x128 Type=UInt16, ColorInterp=Red
Computed Min/Max=664.000,48804.000
Band 2 Block=128x128 Type=UInt16, ColorInterp=Green
Computed Min/Max=3688.000,30697.000
Band 3 Block=128x128 Type=UInt16, ColorInterp=Blue
Computed Min/Max=19.000,38786.000
Band 4 Block=128x128 Type=UInt16, ColorInterp=Undefined
Computed Min/Max=2521.000,39791.000

So I've stretched the contrast with GDAL

gdal_translate test.vrt out.tif -scale 664 48804 0 65535 -exponent 0.5 -co COMPRESS=DEFLATE -co PHOTOMETRIC=RGB

but the result has even less contrast: Link (all export from QGIS).

What went wrong?


texts: Setting bounds for image overlay in Leaflet

texts:
I'm using LeafletJS map. I need to use an image overlay. I have done it like so:

onMapReady(event: any) {
   this.layers.push(imageOverlay(this.imgSrc + "map.png", [[35.334847, 25.328072], [35.337847, 25.339072]])); // This is the overlay
}

I can see the small image overlay as below. But how can I show it the whole screen? I think the problem here is set bounds. But how can I set it? I have center  "lat": "35.335847",
        "lon": "25.338072", If this is the center (i.e. with arrow) then how can I show the image entire page? Any clue about how to set the bounds and find out those lat and lon for that?


texts: Reconciling Differences Between Python List Sort and ArcMap Table Sort

texts:
I have a feature class being used in Data Driven Pages, which sorts as follows:

APPLE LARRY
JOHNSON SHIRLEY
KENT-STATE-GOODS LLC
KENT TIMOTHY
KENT VANESSA
KENT WAYNE
SMITH JEEVES

I also have the Data Driven Pages run via a script, and because we don't yet have ESRI Production Mapping, table generation has to be done as part of this script.  When this script runs, a Python list is created containing these names.  My list parcelOwners[] is then sorted using parcelOwners.sort(), but I get a different sort order amongst the KENT's because of the hyphen.  That sort order is as follows:

APPLE LARRY
JOHNSON SHIRLEY
KENT TIMOTHY
KENT VANESSA
KENT WAYNE
KENT-STATE-GOODS LLC
SMITH JEEVES

As you could imagine, the difference in sort order causes problems when the tables generate for a given set of pages amongst the KENT names.  Page 3 gets the table data intended for Page 4, 4 for 5, 5 for 6, and 6 for 3.  Is there a way to change how Python sorts my list when it encounters hyphens to match how ArcMap sorts the same data?


texts: Google Earth Engine, get a List of Bands

texts:
I want to Automate Band selection for all bands. Therefore a List of the Bands is needed. This didn´t work:

ee.List(image.get('bands'))


texts: Showing roads that are below 10 ft elevations using ArcGIS Desktop?

texts:
I have point elevation data and road layers for a region.

How do I show the road sections that are below 10 ft elevations?

I am using ArcGIS Desktop 10.3.1.


texts: Scaling vector field marker arrow based on attribute?

texts: I have discovered how to create vector lines to indicate flow direction and velocities. With a "vector field marker" symbol layer type, the lines lengths vary based on the data attributes. However, when using an arrow as a symbol layer type, I can only seem to get the line to scale and not the arrow head. I would like to scale the arrow head proportionally with the line based on the data attributes. Is there a way to do this?


texts: Selection/Extraction of features based on location AND attribute

texts:
My data consists of a point layer (iceberg positions containing their acquisition date as attributes) and a polygon layer (areal extends of DEMs also containing the acquisition date).
My aim is to select/extract the features based on their same location (the iceberg positions overlapping with the DEM extends) AND under the condition that they possess the same acquisition date.
Unfortunately I couldn't figure out how I can solve this problem.
What I tried is:

Extracting by location, which only fulfills the condition that their locations overlap
To figure out how to run an extraction/selection based on the attributes of different layers, which I haven't figured out yet


texts:
I'm using the Carto API for a few functions and I'm looking for a way to return geocoded information similar to how Google does it here. With Google, you make a call with the Address like: https://maps.googleapis.com/maps/api/geocode/json?address=1600+Amphitheatre+Parkway,+Mountain+View,+CA&key=YOUR_API_KEY and it returns a JSON like:
{
   "results" : [
      {
         "address_components" : [
            ...
         ],
         "formatted_address" : "1600 Amphitheatre Pkwy, Mountain View, CA 94043, USA",
         "geometry" : {
            "location" : {
               "lat" : 37.4267861,
               "lng" : -122.0806032
            },
            "location_type" : "ROOFTOP",
            "viewport" : {
               "northeast" : {
                  "lat" : 37.4281350802915,
                  "lng" : -122.0792542197085
               },
               "southwest" : {
                  "lat" : 37.4254371197085,
                  "lng" : -122.0819521802915
               }
            }
         },
         "place_id" : "ChIJtYuu0V25j4ARwu5e4wwRYgE",
         "plus_code" : {
            "compound_code" : "CWC8+R3 Mountain View, California, United States",
            "global_code" : "849VCWC8+R3"
         },
         "types" : [ "street_address" ]
      }
   ],
   "status" : "OK"
}

I'd like to do the same with Carto, but I'm only finding how to parse the information through a SQL Update, but not return it on a GET call.


texts:
I'm having trouble adding a WMS to XMap 8, and I'm curious if it can be done at all.

I've been using Global Mapper for geologic research for quite a while, which has the ability to download data from inside the program, from a number of sources (World Imagery is of particular use to me, as well as NE US 10m DEM and NE US Topo).

The company I work for now uses XMap 8, which seems to be a little more robust than Global Mapper, at the cost of the easy data download (or so it seems). I'm in the process of learning Xmap, which must be used due to license restrictions.

Essentially, I need to be able to zoom to a location in XMap, download a topo map, imagery, and elevation and work with the data (the tools are able to do everything I need at this point, but getting the data is a bit troublesome).

If there's a simple workaround or quick way of plopping the maps into XMap, I would be more than happy to hear them.


texts: Interpreting GEE unmixing SMA classification output results

texts:
I would like help interpreting unmixing SMA (Spectral Mixture Analysis) output results in GEE.
I have identified endmembers by class as the following:
4 classes: 0 = forest; 1 = savanna; 2 = shadeWater; 3 = bareground
The following is the code snippet where I merged all the training points together.
//////////////////////////

// merge all the points in the same feature collection
    var points = forest.merge(savanna).merge(shadeWater).merge(bareground)
    print(points)

// extract the bands values using the points defined above
    var endmembers1 = greenestComposite.first().sampleRegions({
  // Get the sample from the polygons FeatureCollection.
  collection: points,
  // Keep this list of properties from the polygons.
  properties: ['class'],
  // Set the scale to get Landsat pixels in the polygons.
  scale: 30
});
print(endmembers1.toList(4), 'endmembers1')

The following is the return of the pixel information:
var unmix = image
  .clip(DAD1geom)
  .select(bands)
  .unmix(endmembers2.getInfo(), true, true);//the trues constrain the proportion results between 0 and 1

  return unmix.addBands(image)
};
////////////////////////////

Here are my results for one pixel after the mapUnmix() algorithm:
FINAL: Image (13 bands)
band_0: 0.19283048751277443
band_1: 0
band_2: 0.2873728238557205
band_3: 0.5197966886315166
B1: 279
B2: 502
B3: 529
B4: 2122
B5: 2028
B6: 1085
B7: 2969
NDVI: 0.6009052991867065
DOY: 16

I am interested in the band_ values because it tells me the proportional values of each class, band for that particular pixel.
What I don't understand is whether each band represents each endmember?
Does band 0 represent forest, band 1 represent savanna, band 2 represent shadewater and band 3 represent bareground for each pixel?
Or do the bands differ based on which endmember I have assigned and am looking at?
I am confused because those bands don't tell me what I expected. Therefore I can't truly tell what the band values mean.
What does the output of bands return?


texts: Complete lookup table for German admin codes in Geonames data? (admin4_code)

texts:
I am working with both the postal and Geonames data for Germany. I noticed that while the Postal data may be used as a lookup for ADMIN_CODES 1-3 (ADM1,ADM2,ADM3), it does not mention ADM4 at all.

Does anyone know where to find the lookup table for 'admin4_code' values?

According to http://forum.geonames.org/gforum/posts/list/538.page, this admin code references a 'Gemeinde', an admin field I would like to include in my work, if possible.

The German dataset I'm working with: http://download.geonames.org/export/dump/DE.zip

The German dataset I'm using as a lookup table for ADMIN CODES: http://download.geonames.org/export/zip/DE.zip

Both zipped folders contain a 'readme.txt' file that outlines the fieldnames in the data.

I am also tracking this question on the Geonames forum: http://forum.geonames.org/gforum/posts/list/2623.page


texts: Importing external bookmark (*.dat) file using ArcGIS Engine?

texts:
Does anyone know how to use ArcObjects (specifically ArcGIS Engine for me) to import a bookmark (saved as a .dat) that was created and saved externally using ArcMap?

I have an ArcGIS Engine application that uses the ControlsMapManageBookmarksCommandClass, but this does not have Load/Save options that the ArcMap bookmark control does.

I've searched all documentation and forums with no luck.


texts: Value Error when adding adding a layer in ArcGIS Pro arcpy script

texts:
In a script, I'm creating a feature layer and trying to add it to a map in ArcGIS Pro 2.3.0. When I try to add it though, I get a ValueError. Any ideas why and how to fix it? m is the map object. The layer is created.

FWIW, I'm running the script in Spyder, Python 3.6.6, which I installed (maybe correctly, maybe not) in the default ArcGIS Pro environment. It seems to recognize the arcpy module, although there's no intellisense.

aprx = arcpy.mp.ArcGISProject(prj)
m = aprx.listMaps()[0]
target = m.listLayers('Cities Villages Twps')[0]
try:
    target.definitionQuery = "NAME = 'Chicago'"
    arcpy.CopyFeatures_management(target, 'GPP')
    arcpy.MakeFeatureLayer_management('GPP', "Area501")
    print(type("Area501"))   # str
    print(type(Area501))     # NameError: name 'Area501' is not defined
    print(type(Area501[0]))  # NameError: name 'Area501' is not defined
    m.addLayer(Area501[0])   # NameError: name 'Area501' is not defined

except Exception:
    print(arcpy.GetMessages())
    print('\n')
    tb.print_exc(limit=1, file=sys.stdout)

del aprx


texts: Any updates on problems with loading ESRI FileGDBs in QGIS 2.2?

texts:
There are many posts in different forums about people having trouble reading File GDBs in QGIS 2.2, but no real solutions. Is there any word on a resolution to this problem?

I am using QGIS 2.2 (installed via the OSGeo4W online installer) on Windows 8.1 64-bit. I can confirm that the gdal-filegdb (version 1.11.0-1) library is installed, as reported by the OSGeo4W setup.

The FileGDB option does not appear at all in the Directory source type list, and if I go and choose the .gdb folder anyway I get an error about it not being a valid data type.

And please don't suggest compiling QGIS from scratch with the gdal-filegdb library linked -- that might fly on Linux but it is not a practical alternative on Windows, even assuming I did know how to do it.


texts: Split line across multiple polygons preserving order

texts:
I have polygons stored in the database as (Multi)Polygons using PostgreSQL (13.1) and PostGIS (3.0.2) with a srid of 4326 (aka states). A user generated line (LineString) (aka route) will be used to find all intersections of those polygons along that line.
I am able to get all of the intersection points using the following query (using this route / linestring):
SELECT states.*
  , (ST_Dump(ST_Intersection(ST_Boundary(states.boundary_polygon::geometry), :route::geometry))).geom AS point
FROM states
WHERE ST_Intersects(:route::geometry, states.boundary_polygon::geometry)

The issue is, I can't seem to find a way to keep the intersection points in order along the route. I have tried to use order by distance from the start point, but the issue is I need both points per each state, thus causing the order to possibly be wrong.
For example, I have a route that splits states A, B, C. I need the start point to border crossing of A. I then need the border crossing of B to C. Then C to the end point.
The reason I need this is to calculate distance, which I currently handle in a separate query as I'm not sure how to combine them into one query yet (I'm a total noob to postgis and my main role isn't a DB admin).
I have tried:
ORDER BY ST_LineLocatePoint(ST_GeomFromText(:route_str, 4326)::geometry, (ST_Dump(ST_Intersection(ST_Boundary(states.boundary_polygon::geometry),ST_GeomFromText(:route_str, 4326)::geometry))).geom)

and
ORDER BY :start_point <-> point
  , :start_point <-> states_with_breakpoints.boundary_polygon::geometry
  , :end_point <-> point
  , :end_point <-> states_with_breakpoints.boundary_polygon::geometry

but still don't get the correct order.
Expected order: start point, A border, B border, C border, end point
Output order: start point, B border, A border, C border, end point
I'm very much a beginner in this and have been searching the webs and using postgis docs to get as far as I have now. Any help would be appreciated.
Env Info: PostgreSQL 13.1, POSTGIS="3.0.2 2fb2a18" [EXTENSION] PGSQL="130" GEOS="3.8.1-CAPI-1.13.3" PROJ="7.1.1" LIBXML="2.9.10" LIBJSON="0.15", [RoR] Ruby 2.5.3, Rails 5.2.4.2 (going to be upgrading to 5.2.4.4 very soon), rgeo gem 2.1.1

Just so you have it, this is the query I'm using for distance:
SELECT ST_Length(
  ST_LineSubstring(
    route,
    LEAST(ST_LineLocatePoint(route, start_point), ST_LineLocatePoint(route, end_point)),
    GREATEST(ST_LineLocatePoint(route, end_point), ST_LineLocatePoint(route, start_point))
  )::geography
) AS distance
FROM (
  SELECT :route AS route
    , :start_point AS start_point
    , :end_point AS end_point
) AS distance_sub_query

Update 1
My end goal is to have a query that splits the route in each state and gets the distance of each route segment in that state. There could be multiple segments in each state. What i'll need is the border points for each state (where the route crosses) and the length of each segment.
I'm open to other ideas or suggestions if the way I'm going about this is wrong or impossible. Would like the query to be under 1 second if possible.

Update 2
Updated query:
SELECT s.id, s.abv, dmp.geom
FROM   states AS s,
       LATERAL ST_Dump(ST_Intersection(ST_Boundary(s.geom), :route::GEOMETRY)) AS dmp
WHERE  ST_Intersects(s.geom, :route::GEOMETRY)
ORDER BY dmp.path

Results:

Order should be: OH -> WV, WV -> PA, PA -> WV, WV -> MD, MD -> PA, PA -> NJ, NJ -> NY

Update 3
SELECT s.id, s.abv, dmp.geom
FROM   states AS s,
       LATERAL ST_Dump(ST_Intersection(s.boundary_polygon::geometry, :route::geometry)) AS dmp
WHERE  ST_Intersects(s.boundary_polygon::geometry, :route::geometry)
ORDER BY dmp.path

This is now not returning the border points and is still out of order.


texts: Join attributes of 2 points to 1 polyline

texts:
I have a set of lines(polyline) that were created as lines between point pairs. Now I need to know, which points the lines connect. Also I want specific column from points table to copy in line table automatically.
Lets say points contain: ID, CODE, ... And in lines I want to copy CODE column from points.
Example line set: ID, CODE1(code column from first point that created the line), CODE2(code column from second point that created the line)

On the image below I tried code that work for points which are associated with 2 lines but when one point is associated with 3 or more lines I get error which can be seen on the same image. Does anyone see some bug in the code?


texts: Wrong projection on WFS layer with OSM baselayer

texts:
I am trying to add a vector layer with a WFS protocol.

var wfs = new OpenLayers.Layer.Vector("WFS", {
                strategies: [new OpenLayers.Strategy.BBOX({resFactor: 1})],
                projection: new OpenLayers.Projection("ESPG:900913"),
                protocol: new OpenLayers.Protocol.WFS({
                    url: geoserverURL,
                    srsName: "EPSG:28992",
                    maxFeatures: 20,
                    featureType: "FEATURE",
                    geometryName: "the_geom",
                    featureNS: "http://tempuri.org/"
                })
            });
            map.addLayer(wfs);

In my database (postgis) the coordinates are stored as EPSG:28992.
My baselayer is OSM so that should be ESPG:900913.
But the features are not projected on the correct spot.

I tried other projections but the features are not projected correctly that way either.

What am I doing wrong here?

Edit:

I tried changing

projection: new OpenLayers.Projection("ESPG:900913")

to

projection: new OpenLayers.Projection("EPSG:4326")

or

projection: new OpenLayers.Projection("EPSG:28992")

But the features are still nowhere near where they are supposed to be.


texts:
've been working on an app for a presentation and we are using ArcGIS Online to geoenrich Canadian Census Dissemination Areas upon the users click. The only catch here is that some of the polygons are causing the request URLs to exceed the 2048 character limit in URLs. We are using the ESRI Proxy.ashx that you can get from GitHub and so i switched from a get request to a post request only the Post request doesn't seem to recognizes the variables being passed in the form data.

This first link is the request in Raw Format as captured by fiddler

http://imgur.com/xlWfh6F

The second link is the same request in WebForm view

http://imgur.com/3fPoHlL

This is the code i'm using to make the request to the proxy server, the proxy adds the token, after it gets the information from our ArcGIS Online account

var studyAreas = {
    /*    areaType: "RingBuffer",
        bufferUnits: "esriKilometers",
        bufferRadii: [2],*/
    geometry: evt.graphic.geometry
}
var sa = JSON.stringify(studyAreas);
var dem = new Demographics;

var cats = dem.categories;
var av = '[';
for (l = 0; l < cats.length; l++) {
    var curCat = dem[cats[l]]
    for (m = 0; m < curCat.length; m++) {
        av += curCat[m].varName;
        av += ',';
    }
}
av = av.substr(0, av.length - 1);
av += ']';
request.post('/proxy/proxy.ashx?http://geoenrich.arcgis.com/arcgis/rest/services/World/GeoenrichmentServer/Geoenrichment/enrich', {
    data: {
        StudyAreas: sa,
        analysisVariables: av,
        forStorage: false,
        f: 'json'
    }
});


texts: QGIS: Building heights and DTM to STL (or other format that CAD can read)

texts:
I have downloaded a DTM and building heights from UK Digimaps (we have institutional access) and would like to:
Export them to STL or other CAD format. Is it possible to do this in QGIS or do I need multiple programs?
I haven't found any way of exporting to STL (dwg or other) yet.

Seems GlobalMapper is an option: http://www.bluemarblegeo.com/products/global-mapper-formats-elevation.php#STL_StereoLithography_Files
Isn't there a way in QGIS ?


texts: Getting bounding boxes for all polygons in GeoJSON feature collection?

texts:
I have a geojson file containing a FeatureCollection, where each feature contains a geometry of type "Polygon". I would like to replace the polygons to rectangular bounding boxes in a simple batch manner. That being said, I could also live with the geometry remaining of type "Polygon", but then of course with just four nodes consisting of the min and max lat/lon coordinates of the former multi-node polygon.

Is there any simple way to do it, favorably some existing algorithm or web service?


texts: Getting displacement in meters from "XYZ ECEF" or "NMEA" data

texts:
I am new to this field. I am working on a GNSS receiver that gives me live location data in NMEA and XYZ ECEF format.
The receiver is mounted on my rover and I wish to use this data to get the rover's displacement in X and Y direction, in meters.
Is that possible ? How should I go about it?


texts: Table Attributes to be merged from Multiple identical polygons in ArcGIS Desktop?

texts:
I have an external Dataset containing Cadastral data. When the lot is shared between different landowners, there are multiple overlapping polygons each storing the landownership and other important attributes separately.

Every overlapping polygon has the same properyID - common feature.

I can extract the table bring it to an SQL engine and create the queries manually; merging all the landowners (L1,L2,L3 etc) and linked it back to a dissolved table by ID.

Can anybody recommend an easier way to do it on ArcGIS Desktop?


texts: Mosaic to Raster Tool in ArcGIS Darkening Grayscale Raster Output

texts:
I stitched roughly 20 raster files together using the "Mosaic to New Raster" tool. Once these raster files were combined, though, the new output raster seemed to "darken"- the fine details of the grayscale images I am using became much less apparent when they were combined in to one file. The image quality is important because I am doing a cost-distance analysis, where the different shades are representative of elevation. For reference, please take a look at the images below.

Is there any likely reason for this change in picture quality? For reference, when using the Mosaic to New Raster tool I used a Pixel Type of 32-bit-Float, with one band, a mosaic operator of "Last", with a TIFF file type. I am also using ArcGIS version 10.1.


texts: Central meridian in Lambert conformal conic distorting distance?

texts:
Does a different central meridian in Lambert conformal conic projection with two standard parallels make a difference in line length?

I have projected a line in GCS WGS84 TO Lambert conformal conic projection,
with the standard parallel 1 = 09 and standard parallel 2 = 35. I have no idea what should be the central meridian.

The actual length is short in Lambert conformal conic projection system.

The actual Route from kathmandu to Delhi is asshown in figure below from skyvector.com and the distance is 460 nm ..............

I have drawn the same route in arcmap with in the following projection

GCS_WGS_1984
WKID: 4326 Authority: EPSG
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
Spheroid: WGS_1984
Semimajor Axis: 6378137.0
Semiminor Axis: 6356752.314245179
Inverse Flattening: 298.257223563

I created the following projected coordinate system according to jeppesen chart.

Asia_Lambert_Conformal_Conic_09and 35
Authority: Custom
Projection: Lambert_Conformal_Conic
False_Easting: 0.0
False_Northing: 0.0
Central_Meridian: 70.0
Standard_Parallel_1: 9.0
Standard_Parallel_2: 35.0
Latitude_Of_Origin: 22.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563

Jeppesen chart with the projection system Highlighted in Red

I projected my line from WGS84 GCS to above Defined projected coordinate system and when I calculated the length of the line in ArcMap the line is 10 NM shorter than the actual one.
I believe with the same projection as Jeppesen Chart, I should be able to get the length of line exactly as Jeppesen chart. I don't know where I am mistaken.


texts:
I am trying to stop the edits using the python console of the QGIS and the code QgsTrackedVectorLayerTools.stopEditing(layer,False) where layer is the layer in the edit mode.

I am receiving the following error:
NameError: name QgsTrackedVectorLayerTools is not defined

Does anyone know the right syntax for the QgsTrackedVectorLayerTools function?


texts: Organizing Data across multiple projections

texts: My department is in the process of migrating data to a more consolidated structured series of GDB's. We are working toward something along the lines of of having 5 different GDB's and placing all of our data in those. For example instead of having multiple one-off GDB's for leasehold we will have it all in one GDB. The issue would be that we have leases in 5 different projection zones (NAD 27 State Plane) thus creating the need to have our GDB in one of those projection zones or having the GDB not have a projection. Then we would just project the MXD's that the data is used in. I am trying to avoid moving a lot of data and then later down the road finding issues with this approach. Are there any issues or a perhaps a better way to approach this issue?


texts: 'Faulty' Landsat-8 Images (Collection 1 Tier 2) for Antarctica?

texts:
Whilst working with the Landsat collections across Antarctica I have noticed that recent L8 images (from late 2017-present) appear to not be 'correct' in GEE. Images from this time period appear black in the viewer, and when you inspect the pixel values they are all negative (no data). However, images before this time period appear to be fine.

I have checked the Landsat repository in Earth Explorer, and the images seem to hold data, so I am wondering if there is a fault within GEE, or whether I have missed a trick here?

Attached is a simple script, showing a 'working' image and a 'faulty' image. Any advice would be great!

https://code.earthengine.google.com/14c0b1c5bf61bc40627a3133cd9c1116


texts: Calling external plugins downloaded via Plugins through a separate out of QGIS Python Script

texts:
I was wondering if I could produce webpages daily by automatically creating a standalone python script, outside of QGIS.
My webpage should ideally contain various layers and an interpolation on top of them.
I was wondering if without opening QGIS at all, I could use the QGIS2Web Plugin, as well as the Openlayers Map Plugin and call the requisite functions from these in a Python Script that I can schedule. I would want all of this to work on a Virtual Machine - so can't really call QGIS.

Note: I did come across various similar questions but there isn't any specific answer to my question.

I have so far been able to pretty much get most of my work working on the Python Console inside of QGIS. (I am using 2.18) - I was able to find the functions for basically reading shapefiles, performing interpolation, and then to call the function for QGIS2Web Plugin. What I aim to do is to essentially automate all of this, with different data on a daily basis.

Edit: I understand I haven't been able to put across the issue. Let me clarify.

I have a shapefile of a particular location, and data for a few points inside this shapefile's extent. I would like to interpolate over this region and produce a heatmap. I realised that using the Openlayer Plugin, I could add a google satellite image map layer to the current layers consisting of the layers for interpolation, data, as well as the initial shapefile, making the heatmap look more realistic as it would now be on top of an actual geographical area. Then I was just using the QGIS2Web tool to make HTML pages out of them.

What I really want to do is to automate this process, since I will be receiving the data real-time through sensors.
Till now I have been successful enough in running successfully the processes of reading the data, the shapefile and performing interpolation, along with using QGIS2Web to create the HTML, through the in-built Python Console present in QGIS.

I was wondering if I could just do all of the above processes by using a Python Script outside of QGIS, by just appending the corresponding directories to which these belong.


texts:
I am fairly new to spatial analysis and I am trying to create a map of the spatial distribution of pollution in the UK using publicly available data retrieved from https://uk-air.defra.gov.uk/data/pcm-data.
The files are all downloadable as .csv and provide the x and y coordinates for the centroid of each grid cell.
I tried to import them on QGIS 3.10 in the following way: Layer > Add Layer > Add Delimited Text Layer... Once I select the file, QGIS automatically recognizes the centroids as my x and y coordinates. However, when I try to plot them, my pollution concentrations do not overlap with UK geographical boundaries.
My ultimate goal is to join attributes from different datasets by location so I would need help understanding how to properly import my pollution csv files so that they overlap with my other shapefiles.


texts: Redefine (or reconfirm) Node and Vertex in 3D-GIS

texts:
Sorry if this sounds too broad and opinion-soliciting...

I think the difference between Node and Vertex had been discussed before (for instance) Difference between node and vertex in GIS?, good and detailed discussions between renowned GIS leaders are open for us (me) to learn.

And my take is (sorry again if this is too simplistic... and I can be totally wrong...) "Vertex" is for location (xy) and "Node" includes its topology.

Is it still valid, given we now have 3D-GIS more closer than before? I mean; Vertex is the XYZ locator?


texts: Which SAGA versions does QGIS 2.18.10 support?

texts:
Beginning with QGIS 2.18.10 the SagaAlgorithmProvider do not longer have a dictionary called supportedVersions. Previously, it looked like this:

supportedVersions = {"2.1.2": ("2.1.2", SagaAlgorithm212),
                     "2.1.3": ("2.1.3", SagaAlgorithm213),
                     "2.1.4": ("2.1.4", SagaAlgorithm214),
                     "2.2.0": ("2.2.0", SagaAlgorithm214),
                     "2.2.1": ("2.2.0", SagaAlgorithm214),
                     "2.2.2": ("2.2.2", SagaAlgorithm214),
                     "2.2.3": ("2.2.3", SagaAlgorithm214),
                     "2.3.0": ("2.3.0", SagaAlgorithm230),
                     "2.3.1": ("2.3.0", SagaAlgorithm230)}

Now it simply says:

if not version.startswith('2.3.'):
  ProcessingLog.addToLog(ProcessingLog.LOG_ERROR,
                         self.tr('Problem with SAGA installation
                                  unsupported SAGA version found.'))

So I would assume that now only SAGA versions starting with 2.3 are supported. However, in my QGIS GUI the Log Messages Panel tells me that SAGA version 2.3.1 is not supported (interestingly, you can use SAGA 2.3.1 in a custom application like RQGIS). So which SAGA versions are supported in QGIS 2.18.10?


texts: Seeking alternatives to Geodatabase to describe topological relations?

texts:
According to switchfromshapefile.org, there is no way to describe topological relations in a shapefile. I am aware that this can be done using a Geodatabase and I have heard of TopoJSON.

If there are any other alternatives which GIS Software works with them?


texts: Bounding box in Dijkstra algorithm (many to many)

texts:
I apologize for the delicate duplication of my TOPIC but I am still standing.

I have the table, which has information about pairs between source places and points (without geometry).
I have used this code so far:

drop table if exists dijkstra;
create table dijkstra as (
SELECT a."id_pk"
, a.mz_id
, a.poi_id
, case when a.nearest_mz!=a.nearest_poi then SUM(a.cost)
else 0
end AS "tmp_cost"
FROM
( SELECT "id_pk", mz_id, poi_id, nearest_mz, nearest_poi
, (pgr_dijkstra( 'select id::integer, source, target
, ST_Length(geom) as cost from ways',
 "nearest_mz", "nearest_poi", false )).*
FROM pairs_mz_poi_node where id_pk <100 -- limit for test )
AS a
GROUP BY a."id_pk"
, a.mz_id
, a.poi_id
, a.nearest_mz
, a.nearest_poi)
;

It works fine, but now I have too many pairs and I have to calculate it faster than before.

I'm trying to add bounding box in my code.

I used this codes:

drop table if exists test;
create table test as (
SELECT a."id_pk"
, a.mz_id
, a.poi_id
, case when a.nearest_mz!=a.nearest_poi then SUM(a.cost)
else 0
end AS "tmp_cost"
FROM
( SELECT "id_pk", mz_id, poi_id, nearest_mz, nearest_poi
, pgr_dijkstra ( 'SELECT gid AS id, source, target,
ST_length(geom) AS cost FROM ways
 WHERE geom && ST_Expand(
(SELECT ST_Collect(the_geom) FROM ways_vertices_pgr d
join pairs_mz_poi_node b
on b.nearest_mz=d.id
join pairs_mz_poi_node c
on c.nearest_poi=d.id
WHERE id IN( b.nearest_mz,c.nearest_poi ) )
             nearest_mz i nearest_poi
, 1)  ',
          "nearest_mz", "nearest_poi", false )).*
FROM pairs_mz_poi_node where id_pk <100 )
AS a
GROUP BY a."id_pk", a.mz_id
, a.poi_id
, a.nearest_mz, a.nearest_poi);

And this:

drop table if exists test2;
Create table test2 as (
SELECT a."id_pk"
, a.mz_id
, a.poi_id
, case when a.nearest_mz!=a.nearest_poi then SUM(a.cost)
    else 0
    end AS "tmp_cost"
FROM
( SELECT "id_pk", mz_id, poi_id, nearest_mz, nearest_poi
, (pgr_dijkstra( 'WITH
buffer AS (SELECT id, ST_Buffer(the_geom, 1500) AS geom FROM ways_vertices_pgr)
                select d.id::integer, d.source, d.target
                , ST_Length(d.geom) as cost from ways d join buffer
                on buffer.id=d.source or buffer.id=d.target
where d.geom=buffer.geom'
                , "nearest_mz", "nearest_poi", false )).*
FROM pairs_mz_poi_node where id_pk <100 )
AS a
GROUP BY a."id_pk", a.mz_id
, a.poi_id
, a.nearest_mz, a.nearest_poi
);

I found this topic but I don't know, how I can change ID number to nearest_mz and nearest_poi from my table.


texts: QGIS Semi-automatic classification plugin install error on OS X 10.10: matplotlib.backends.backend_qt4agg

texts:
I have tried repeatedly to install the QGIS Semi-automatic classification plugin via both (1) Plugins > Install & Manage Plugins and (2) manually based on other comments I found here, but I can't get it to work. I continue to get the error:

The plugin is broken. Python said:
No module named matplotlib.backends.backend_qt4agg

I have checked and have installed matplotlib, GDAL, scipy, numpy via macports (QGIS was installed via macports as well, but I have uninstalled, reinstalled, and rebooted a couple of times at this point).

I am not sure how else to proceed.


texts:
How is the st_length function stored in an Oracle geodatabase? Is it in a:

PL/SQL function
Separate .dll library
Somewhere else?

I've noticed that the sde.st_length function works, even when the st_shapelib.dll library is not configured. I can call sde.st_length without issue: SELECT SDE.ST_LENGTH(SHAPE) FROM ROADS_SAMPLE. So clearly, st_length is not stored in the st_shapelib.dll library.

Knowing this, my next guess would have been that the st_Length function would be stored in PL/SQL in the database somewhere. However, when I look in the Functions  folder in the database (using SQL Developer), I don't see any sign  of the st_length function.

Update:

I've found two functions in the database that pertain to st_length:

ST_LENGTH_F:

Located: Packages.ST_GEOMETRY_OPERATORS Body.st_length_f

Function st_length_f(prim SDE.st_geometry)
Return number
IS
Begin

  If prim IS NULL OR (prim.numpts = 0 and prim.len = 0) THEN
    Return NULL;
  End If;

  Return(prim.len);

End st_length_f;

ST_LENGTH:

Located: Types.ST_GEOMETRY.ST_GEOMETRY Body.st_length

member Function st_length Return number IS
  Begin
    Return len;
   End st_length;

But I'm not clear on which of these two functions is the one in question. The more obvious choice is st_length_f, but it's got a different name than st_length, so I can't imagine that it is the function that is being called.


texts:
How do I edit metadata in QGIS?

I found the plugin Metatools which is experimental and did not work for me.

In Metatools I always get the information

  ...unsupported metadata standard! Only ISO19115 supported now!

I then exported the metadata of the shapefile to xml with the same name as the shapefile, expecting to get the metadata in ISO standard but that is still not working.

My System:
Linux Mint 17; QGIS 2.6.0 Brighton; Metatool 0.3.0


texts: Examples of Google styled Maps? (including JSON)

texts:
After watching a great video on creating different map styles with Google Maps, I have been having some fun playing around with the Google Maps Style Wizard.

I was hoping to find a gallery of user submitted Google Maps Styles, with accompanying JSON code so that I could pinch their hard work & cartography skills, and use on my own applications.

Has anyone come across either a collection of styles (+JSON) or at least some great examples of google maps styles, where the JSON is also available?

The closest I could find was this.


texts: Getting coordinates of pixels inside QgsRasterBlock

texts:
I am trying to implement a function that filters some points around a focus point in QGis, based on a raster map, and then does some distance calculations with the result. For this, I need to convert between geo-coordinates and raster layer point indices.

So, let's say raster is a QgsRasterLayerObject. I know I can get the data array like

provider = raster.dataProvider()

raster_extent = provider.extent()
raster_width = provider.xSize()
raster_height = provider.ySize()
no_data_value = provider.srcNoDataValue(1)

block = provider.block(1, raster_extent, raster_width,
            raster_height)
if block.isValid():
  for x in range(raster_width):
    for y in range(raster_height):
      value = block.value(x, y)
      if value != no_data_value:
        ...

I have read the code of QgsRasterDataProvider.identify and it seems to me that there is only very simple arithmetic computation deriving the array index from the coordinates, which in Python looks like this.

extent = provider.extent()
xres = (extent.width()) / provider.xSize()
yres = (extent.height()) / provider.ySize()
col = int((location.x() - extent.xMinimum()) / xres)
row = int((location.y() - extent.yMinimum()) / yres)
value = block.value(col, row)
if value != no_data_value:
  ...

Are there any data providers of rectangular data for which this code would lead to wrong, but not obviously apparently wrong results, or can I use this and its obvious inverse in my function?


texts: how to split a line by a point which is on the line but not in the linestring itself?

texts:
i need some help with the following:

i have a linestring - geom
i have $lat,$lon (using some php syntax)

i have an osm table with the roads geom, i am trying to get closest point, and then split the road into two parts, by that closest point.

after using the lat,lon and st_closestpoint to retrieve the closest point on that line.
the point which will be returned is not necessarily exist in the original linestring points.

ST_ClosestPoint (geom_way, ST_SetSRID (ST_MakePoint($lon, $lat), 4326)  ) as closest_point

after that i would like to split the geom(linestring) into two lines by closest_point,
and this is where i need help.

i tried using

 st_asgeojson (ST_Split (geom_way, ST_ClosestPoint (geom_way, ST_SetSRID (ST_MakePoint($lon, $lat), 4326)  )))

but i believe that because the point i use is not on the geom, its not spliting and always returns the entire geom.

if i use a point which exist in the line string, it works well.(but not what i need)

how can i solve this? any solution even if with totally different approach will do.

thanks for your help!


texts: Can't find the ROI creation dock in SCP Toolbar

texts:
I'm following the tutorial http://semiautomaticclassificationmanual-v4.readthedocs.io/en/latest/Tutorials.html using QGIS 2.16.3 on Ubuntu 16.04. In the section 3.1.4, a shapefile is created using ROI Creation Docker. But I can't locate this in my SCP Toolbar.

The closest I got was the ROI Creation section located under Classification Dock, but even it doesn't allow you to add shapefiles.

Is this because I'm using an incompatible Plugin version or has this capability been replaced by some other feature?


texts: STIntersects and STContains fails on parameter

texts:
I am new to GIS and using GIS in SQLServer, so apologies if this is a rudimentary question, but I cannot find the answer so far.  I am trying to use STIntersects() or STContains to determine what state a given lat/long combination is in.  My table named SpatialTable has polygons for all the states and my example below is coded for a point inside my polygon named AL02.

I need to understand why this code (cloned from example code on the Microsoft TechNet pages) works fine:

DECLARE @g geography;
DECLARE @h geography;
select @g=GeogCol2 from SpatialTable where name = 'AL02';
SET @h = geography::Parse('POINT(-88.09 31.24)');
SELECT @g.STIntersects(@h);

and this code does not:

-- this will not work as it sees my parameter as nvarchar(max)
--error message is: Cannot call methods on nvarchar(max).
SELECT a.name
FROM SpatialTable a
WHERE a.GeogCol2.STIntersects(geography::Parse('POINT(-88.09 31.24)')) = 1
and a.name = 'AL02';

I also tried STGeomFromText instead of Parse, but still get the nvarchar error.

SELECT a.name
FROM SpatialTable a
WHERE a.GeogCol2.STIntersects(geography::STGeomFromText('POINT(-88.09 31.24)',4326)) = 1
and a.name = 'AL02';

Even if I pass @h into the SQL query with @h declared as geography, the error shows that the database is treating the parameter as nvarchar.  How can I dynamically load my lat and long without passing nvarchar?

Any help would be greatly appreciated.  Thanks.


texts: Positional changes bewteen MapInfo 7 and 10

texts:
I have just come across something bizarre and I can't explain it. I have programs for creating grid squares and for converting points to polygons based on an OS grid reference [GB BNG]. Recently we upgraded to version 10 having previously been on version 7. Some data needed updating based on a new protocol for the disclosure of a protected species. The data, using a new easting and northing was given a point and then run through the same polygon generating routine which creates a polygon based on the accuracy of the grid reference.

When projected on current data and grids the newly created data is not aligned correctly. On further investigation all our data appears west 5 meters and 222 metres south of where it should be in MapInfo 10 but is in the correct position in MapInfo 7.

I cannot see what it is that could be different or changed. All the data created in MapInfo 7 appears in the correct position but if it is generated in MapInfo 10 it is not. We've checked alignment against on-line OS data and what we seen in MapInfo 7 is correct!

Some have suggested that MapInfo 10 is more accurate but a metre is still a metre and 10 metres from the projection origin is still the same 10 metres.

Any thoughts?


texts: Struggling with enabling a PostGIS database on Mac

texts:
I have installed PostGIS on the mac using the binaries available on http://www.kyngchaos.com/software:postgres

Postgres v 9.1.3

PostGIS v 2.0

After that I want to enable postgis for my database. I tried using CREATE EXTENSION query. But it error'ed out saying that the extension is not available. I then tried

sudo -u postgres psql -d [your database] -f /usr/local/pgsql/share/contrib/postgis-2.0/postgis.sql
Password:
SET
BEGIN
psql:/usr/local/pgsql/share/contrib/postgis-2.0/postgis.sql:69: ERROR:  could not access file "$libdir/postgis-2.0": No such file or directory

How can i make $libdir expand into the path which i found by executing

$ pg_config --pkglibdir
/usr/lib/postgresql

The /usr/lib/postgresql does NOT contain postgis-2.0.

Digging around for a few more hours I found a hacky solution:

export PATH=/usr/local/pgsql/bin:$PATH

This makes the pg_config magically change to:

$ pg_config --pkglibdir
/usr/local/pgsql-9.1/lib

And then I hardcode that path into the postgis.sql file. But I want to know why $libdir doesnt resolve automatically.


texts: Using v.distance python script parameter with PyQGIS?

texts:
This question is an extension of Copying code generated when you put data into GUI via QGIS Plugins?

I have a Python file I want to run on two loaded layers in QGIS. The layers are TestPhotos (point) and TestObj (point). The tool I want the Python file to run is v.distance.

The column in TestObjects I use in v.distance is OBJ_ID

The Python file is called Test1.py. I have this code:

import processing

processing.runalg('grass7:v.distance', 'TestPhotos','point,line,area','TestObj','point,line,area',-1.0,-1.0,'to_attr','OBJ_ID','OBJ_ID',None,-1.0,0.0001,'nearestobj','distance')

I get the following output using the QGIS Python console:

execfile(u'C:/OSGeo4W64/Test1.py'.encode('mbcs'))
Error: Wrong parameter value: None

I believe it is associated with the v.distance GRASS_REGION_PARAMETER <ParameterExtent> parameter. Any idea how to deal with this to get the Python script working?

EDIT1 - responding to comment

When I typed in processing.alghelp('grass7:v.distance'), i got this:

>>>import processing
>>>processing.alghelp('grass7:v.distance')

    ALGORITHM: v.distance - Finds the nearest element in vector map 'to' for elements in vector map 'from'.
        from <ParameterVector>
        from_type <ParameterString>
        to <ParameterVector>
        to_type <ParameterString>
        dmax <ParameterNumber>
        dmin <ParameterNumber>
        upload <ParameterString>
        column <ParameterString>
        to_column <parameters from to>
        GRASS_REGION_PARAMETER <ParameterExtent>
        GRASS_SNAP_TOLERANCE_PARAMETER <ParameterNumber>
        GRASS_MIN_AREA_PARAMETER <ParameterNumber>
        from_output <OutputVector>
        output <OutputVector>

EDIT2

When I try 'Test.py' with the following code:

import processing
TestPhotos = QgsVectorLayer('c:\OSGeo4w64\TestFiles\TestPhotos.shp', 'TestPhotos', 'ogr')
TestObj = QgsVectorLayer('c:\OSGeo4w64\TestFiles\TestObj.shp', 'TestObj', 'ogr')

extent = TestPhotos.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

processing.runalg(\
'grass7:v.distance', \
TestPhotos,\
'point,line,area',\
TestObj,\
'point,line,area',\
-1.0,-1.0,\
'to_attr',\
'OBJ_ID',\
'OBJ_ID',\
"%f,%f,%f,%f" %(xmin, xmax, ymin, ymax),\
-1.0,0.0001,\
'C:\OSGeo4W64\TestFiles\Outputfile0.shp',\
'C:\OSGeo4W64\TestFiles\Outputfile1.shp')

The output files are correctly generated with the above code. But they are not loaded in QGIS.


texts: Auto populate field with text using field calculator if/then

texts:
I've extracted values from a raster to a point file, but now I need a field with the text information from the raster attribute table that corresponds with the value. I've been trying to use the field calculator to do an if/then type function, but I can't get it to work;
Basically, where it says 1201 in the RASTERVALU field I want it to say "Developed, Open Space" in the CN_LEVEL3 Field. I have absolutely no experience with Python.


texts: Debugging DLL (ArcGIS tool) with Visual Studio 2008 Express?

texts:
As per Esri's recommendation I began the process of converting my existing VBA customization into VB .Net Express edition. I was surprised to realize that debugging a DLL is not an option in the express edition.

"If you are using the Visual C# Express Edition or Visual Basic Express Edition, you cannot debug DLLs because these versions do not have the option to start an external application for debugging."

Source: http://help.arcgis.com/en/sdk/10.0/a...000001w2000000

Any workaround other than purchasing VS Pro?


texts:
I have a .shp file from a grid. I read it and calculate the centroids. This data.frame has the latitude, longitude and some numbers in order making reference to the id of the grid's square. Something like that:
df-centroids

With this, I try to transform it to a raster where the values are the id (like they are different numbers I think it will work), but when I plot the raster it dont use the "id" like values, and in the data of the raster object I think the "id" field has dissapeared:
rasterdf1
 and the image after ploting that shows this: plot(rasterdf1)

I dont know where I am wrong or where I am loosing the value field. Can anyone help me, please?

Here I leave my code:

    shapeData1 <- readOGR("C:/Users/b.vcm/Documents/rural_planner_dashboard_v2",'Grid_40Km')
    #The .shp I use
    trueCentroids <- rgeos::gCentroid(shapeData1,byid=TRUE, id=shapeData1$id)
    df_centroids <- data.frame(trueCentroids, id = shapeData1$id)
    spg2 <- df_centroids
    coordinates(spg2) <- ~ x + y
    spg2 <- SpatialPixelsDataFrame(spg2, tolerance = 0.286937, spg2@data) #tolerance is because coordinate intervals are not constant due to degrees distance
    rasterdf1 <- raster(spg2[,'id'])
    plot(rasterdf1)


texts: Unknown gridlines on old Greenland base map!

texts:
Can anyone help me out please?

I'm using an old 1:10,000 base map of Greenland for a mapping project with nice 10cm gridlines already drawn on. However it uses a coordinate system which is unfamiliar to me, it has values of 92 up to 112 on y axis and 124 to 138 on x axis, which increase by increments of 2. My mapping area is around 61 N, 45.15 W, can anyone please tell me what the grid system my map uses represents?

Thanks
Eleni


texts: How can I adjust the zoom level and recenter in QGIS automatically?

texts:
I am using QGIS 3.14 and the Python console (pyQGIS). I will be tracking a line from GPS co-ordinates and I want to be able to adjust the screen and zoom level accordingly. Here is my code:
canvas=iface.mapCanvas() #get a reference of the canvas where maps are displayed

#Drawing Polyline
polyline = QgsRubberBand(canvas, False)  # False = not a polygon
points =[]
points.append(QgsPoint(-1232443,1318350))
points.append(QgsPoint(-1332443,1418350))
points.append(QgsPoint(-1432443,1518350))

polyline.setToGeometry(QgsGeometry.fromPolyline(points), None)
polyline.setColor(QColor(0, 128, 255))
polyline.setWidth(3)

#This is how I am adding points and plotting them. (Just for testing)
points.append(QgsPoint(-1532443,1618350))
polyline.setToGeometry(QgsGeometry.fromPolyline(points), None)
points.append(QgsPoint(-1632443,1718350))
polyline.setToGeometry(QgsGeometry.fromPolyline(points), None)

How can I adjust the screen size and zoomlevel accordingly?


texts: Selecting all polygons along northern edge of shapefile in R

texts:
I have a shapefile of Mexico composed of Admin 2 polygons.

I am trying to create a buffer along the Mexico-US border, in R considering this is a shapefile of Mexico only and I want to retain the admin 2 boundaries after the buffer is created.

So I am trying to find a way of selecting all the bordering admin 2 polygons. The admin 2 polygons are of varying size with variable numbers of neighbours.

@Spacedman provided feedback -- I found this resource, which was helpful, and my code now looks as follows:

#subset to mexico and us only
mexus <- adm0all[adm0all$ADM0_NAME %in% c('Mexico', 'United States of
America'),]

#identify shared border
sf_bord <- mexus %>%
            ms_innerlines() %>%
            as_tibble() %>%
            st_as_sf()

#plot & see if worked
ggplot()+
  geom_sf(data=mexus) +
  geom_sf(data=sf_bord, color = "red", size = 1) +
  coord_sf(datum = NA) +
  theme_void()

This results in the border being isolated(see image above), however when I calculate the border using the following code, my buffer is fragmented/non-continuous, see code and image below:

#transform to project CRS to get units in meters
sf_bord_utm <- st_transform(sf_bord, 6368)

#calculate 5km buffer
mexus_border <- st_buffer(sf_bord_utm, dist = 5000)

Does anyone understand why this is not a continuous buffer?


texts:
I created a new environment in Conda and installed Python 3.6.9 along with certain modules including ArcPy in that new environment. Now, when the script in Spyder tries to import arcpy, I get the following error:
File "~\Local\Programs\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 14, in <module>
    import arcgisscripting

ImportError: DLL load failed: The specified procedure could not be found.

Please note that I don't use ArcGIS Pro, I only installed the latest ArcPy so that I can use it in a particular Conda environment. I use ArcGIS, but I don't use its ArcPy because it's limited to Python 2.7. How can I resolve this issue?


texts:
I'm trying to build .vrt from georeferenced DEM files from Photoscan.

gdalbuildvrt global_dem.vrt dem*.tif

It was working on Ubuntu 12.04(GDAL 1.10.0, released 2013/04/24), but now it fails on Ubuntu 14.04 (GDAL 2.1.0, released 2016/04/25)

ERROR 1: Failed to write .vrt file in FlushCache().

How to fix it?


texts: Inverting raster changes my min and max values

texts:
I am trying to invert a raster in ArcGIS (10.5), and have used the following expression (from esri support page) in the raster calculator:

invert = (("raster_file" - maximum) * -1) + minimum

However, the high and low values of the new raster are slightly different from the original raster:

original low:0.1551740020513535
original high:0.8855689764022827

new low:0.1852769702673
new high:0.84637397527695

Why is this happening, and what can be done?


texts: Creating Hydro Hierarchy in ArcGIS Desktop?

texts:
We recently ran at work into this model created by ESRI for the rivers in the US:

https://maps.esri.com/rc/radial/index.html
Code: https://github.com/richiecarmichael/Esri-Hydro-Hierarchy

In order to make QA checks we are trying to embed this sort of model into a geometric network of sewage which has about 5000 line features.

We tried to use the add-in attached to the code (ExportHydroNetwork) in order to apply it on the network, but it didn't even appear on the add-in manager after connecting to the folder (arcmap 10.3)

What are the parameters that this add-in requires for work?

Does it, for example, need a field that classifies level of each line segment?

Or maybe it requires modification in the arcmap environment?


texts: Unexpected elevation values from SRTM using Zonal statistics

texts:
In QGIS 2.18.18 on Windows 10, I have downloaded the NASA file srtm_05_01.tif that is available at http://dwtkns.com/srtm and the airports.shp file available at http://qgis.org/downloads/data/qgis_sample_data.zip. I am using them in an exercise to calculate zonal statistics.

I follow these steps to create my analysis area buffer layer and then calculate the zonal statistics

Add the airports.shp layer and the srtm_05_01.tif layer to the project
Vector/Geoprocessing/Tools/Variable Distance Buffer
Accept the defaults and run
Save the resulting buffer layer as a shapefile with the CRS set to EPSG:4326 WGS84, the same as the raster file and project
With the zonal_airport layer created, run Raster/Zonal_Statistics
Select srtm_05_01 as raster layer and zonal_airport as polygon layer, elev as the prefix
Select mean, min, max as the stats
Run.

When I look at the added attributes in zonal_airport, I get the following (only 4 of the 76 Alaska airports are in the srtm_05_01 sector)

Except for Port Heiden, it first appeared that the elevmax in the results can be mutiplied by 3 to get ELEV original attribute (perhaps unit of measure diffence: feet and meters or yards). But with Port Heiden it is apparent that I may have missed something.

What am I doing wrong?


texts: How to get the linestring of a way in postgis with pg_snapshot schema

texts:
I have loaded some data from openstreetmap using osmosis and its pg_snapshot schema. This results in a table wayswhich contains a list of nodesreferenced by their ids:

select w.nodes from ways as w where w.tags->'railway'='rail' limit 1;
         nodes
-----------------------
{358482292,358482290}

and in a separate table you have a geometry column for each node:

select st_astext(nodes.geom) from nodes where id = 358482292;
     st_astext
-----------------------------
POINT(6.1253475 45.9259597)

My question what query could give me a LINESTRING for each way created from the node list?


texts: Offline Base Layer using stored KML or GEOJSON file

texts:
Using OpenLayers 2, I want to have an offline baselayer, that is loaded from a KML or GEOJSON file.  I am successfully able to do this if I first load a baselayer from google maps.  For example, this works fine and am able to properly switch between the two layers:

var googlePhys = new OpenLayers.Layer.Google(
      "Google Physical",
      {type: google.maps.MapTypeId.TERRAIN});
var worldGEOJSON = new OpenLayers.Layer.Vector("GeoJSON", {
projection: "EPSG:3857",
strategies: [new OpenLayers.Strategy.Fixed()],
protocol: new OpenLayers.Protocol.HTTP({
    url: "maps/test.geojson",
    format: new OpenLayers.Format.GeoJSON()
}),
isBaseLayer:true,
visibility:true
});

map.addLayers([googlePhys, worldGEOJSON]);

However, changing just one line, this does not work:

map.addLayers([worldGEOJSON, googlePhys]);

I can switch back and forth between the two and the the google layer works fine, but the geojson layer never displays.  Just a blank white screen.  I have tried changing the style of the geojson layer and that works fine in the first example, but not the second.  Any thoughts?


texts: Randomly placing polygon within delimited area in ArcGIS for Desktop?

texts:
I need to compare the number of animal species in a polygon of my study site against the number of species in a polygon of the same shape and size of my study site but located randomly within a specific region.

Does someone have an idea on what tool/procedure could I use in ArcGIS 10.3.1 (full license) to do so?

Basically I just need to randomly place my study site polygon within the area of a larger region delimited by another shapefile.


texts:
I have a functional SQL query which works in Oracle Spatial, but does not work in PostGis.
What works in Oracle spatial:

"CREATE TABLE mestske_komunikace AS
SELECT r.*
FROM lt_cze_5_0_points_sjtsk r, moravskoslezsky_kraj u
WHERE
(SDO_GEOM.SDO_INTERSECTION(r.geom,u.geom,0.05) IS NOT NULL)
AND
(r.urban='1');"

This returns me picture number 1.

I had to change SDO_Interesection to ST_Intersection in Postgis.
But the script below does not works as it should. It shows me Metropolitan Road even in another region. Do you know why intersection does not work in the Postgis or how to rewrite for it to work?
And another question is, it shows me roads in points and i would like to have the roads in lines, is it possible?

CREATE TABLE mestske_komunikace AS
SELECT r.*
FROM lt_cze_5_0_points_sjtsk r, moravskoslezsky_kraj u
WHERE
(ST_Intersection(r.geom, u.geom) IS NOT NULL)
AND
(r.urban='1');

Picture number 2


texts: Merging two adjacent polygons which borders are not touching each other?

texts:
I want to merge adjacent single polygons which are apart from each other.

Important note: I cannot use ArcGIS, I need to do this programatically, in any language (R, Python, Ruby, PostGIS) doesn't really matter which it is, it can be a combination of a scripting or compiled code and/or command line tools.

So far I have tried using Turf.js methods convex (hull) and union. Union doesn't grow/shrink the edges to to complete with the neighbouring polygon, but it rather just creates a a geometry collection of all the polygons processed as a single feature. The convex hull creates a single polygon but ignores the inner points, just the farthest and gives a wrong final shape. I've tried using the Eliminate Sliver Polygons... in the QGIS Vector>Geoprocessing Tools menu but the operation fails with a Python error: **

  global name 'QErrorMessage' is not defined
  Traceback (most recent call last):
    File "/usr/local/Cellar/qgis-214/2.14.1/QGIS.app/Contents/Resources/python/plugins/fTools/tools/doEliminate.py", line 96, in accept
      self.eliminate(inLayer, boundary, self.progressBar, outFileName)
    File "/usr/local/Cellar/qgis-214/2.14.1/QGIS.app/Contents/Resources/python/plugins/fTools/tools/doEliminate.py", line 236, in eliminate
      QErrorMessage(self).showMessage(
  NameError: global name 'QErrorMessage' is not defined

3 Adjacent polygons that are suppose to be merged into a single one:

Desired resulting single polygon shape
The expected border is drawn with a red contour.


texts: Using z-index with infowindow and CartoDB CSS?

texts:
I currently have a map where my map div covers 60% of the page and 40% is a fixed  Infowindow. My problem happens with cartodb css and its z-index of div.cartodb-infowindow.
Currently my information will not show in the fix div, however the html elements are there when I click on a point and if I manually go into the debugger and un-check the cartodb-infowindow, the information will show up.

So in my styles sheet I added .cartodb-infowindow and tried setting the z-index to a higher number than my side panel but no luck. I have tried a lower number and even setting the z-index to none. However I cannot override the cartodb css on z-index. (I have also tried using !important as well).


texts: Set QGIS 2.18 as default on Win with both 2.18 and 3.0

texts: Is there a way to have both QGIS 3.0 and 2.18 installed on one computer but have 2.18 as the default Package when opening workspaces. At the moment 3.0 opens all my workspaces which is rather annoying as most of the Plugins and providers I'm using are not available yet in 3.0. It was more of a test download to familiarise myself with the newest version until it is brought up to LTR standard. OS is Windows 10


texts: Using GDAL to ID a pair of GeoTiff Coordinates

texts:
I've got a large library of images with no spatial reference. When I view the GeoTiff properties, I can view a lat/long under a "GPS" heading, but there's no indication of a coordinate system/projection.

My goal here is to identify a coordinate system and some coordinates in the image so that it can be georectified.

I've tried using GDALinfo to pull data, but it's giving everything to me in terms of pixels, not coordinates. GDALinfo also isn't outputting any coordinate system information.

Any ideas?

Also tried this code I found online and found "none".

from osgeo import gdal, osr
ds = gdal.Open(r'SOMERASTER.TIF')
prj = ds.GetProjection()
print prj

srs = osr.SpatialReference(wkt=prj)
if srs.IsProjected:
    print srs.GetAttrValue('projcs')
    print srs.GetAttrValue('geogcs')

https://www.dropbox.com/s/3x0jjq7ytu964c6/IMG_0019_1.tif?dl=0


texts:
I want to plot kind of isarhytmoc points map like in this tutorial, but with population data.
However, unfortunately, I have no enough skill and only came with such code:

require(sp)
require(rgdal)
require(RColorBrewer)
require(ggplot2)
library(plyr)
library(maptools)
library(rgeos)
require(gpclib)
gpclibPermitStatus()
gpclibPermit()

rus<-url("http://www.filefactory.com/file/4h1hb5c1cw7r/n/RUS_adm1_RData")
print(load(rus))
proj4.str <- CRS("+init=epsg:3413 +lon_0=105")
gadm.prj <- spTransform(gadm, proj4.str)

popul <- read.csv2(file="C:\\unempl11.txt", header = TRUE,
        sep = ";",quote = "", dec=",", stringsAsFactors=F)

df <- fortify(gadm.prj, region = "ID_1")
df2 <- merge(df, gadm.prj, by.x="id", by.y="ID_1")

p <- ggplot(df2, aes(x = long, y = lat, group=group)) +
       geom_point(data    = popul,
                  mapping = aes(x=lon, y=lat, colour=abs),
                  size    = 3,
                  alpha   = 0.8) +
       scale_colour_gradient2(name  = "Population",
                              low   = "darkred",
                              mid   = "white",
                              high  = "blue",
                              guide = "colorbar") +
       ggtitle("Population in Russia)")
p + geom_path(data    = df2,
               mapping = aes(x=long, y=lat, group=group),
              size    = 0.125)

Here is the data I use, if anyone can help me with some ideas i will be grateful.


texts: Executing published geoprocessing service for customized GP tool from Python scripting?

texts:
I have published the geoprocessing services for a GP tool.

How can I use this link in python scripting to execute the GP tool?

http://VM01-AG01-NJ:6080/arcgis/rest/services/FeatureClassToGeodatabase/GPServer


texts: How to easily change/edit existing field datatype in shapefiles?

texts:
I am editing attribute tables of multiple shapefiles for a future project. I need the field data types to be very specific and uniform among all the shapefiles. I am running into the problem that the majority of these shapefiles's attribute tables are all different in terms of field name and field data type.

I have been adding a new field, then using field calculator to copy the data from the old field column into the new one. However, one field (StationID) has been giving me trouble because copying station id with letters and numbers (example 113a) will not copy into a float in the new added field. And that is fine, I can manually fix that and include the information in a separate Notes column, but not when there over 1,000 entries in which this occurs.

Is there any command, or tool or function, that could "easily" change the field data type from string/text to float?


texts: Float error in converting cub image from ISIS3 toolbox to Gtiff using GDAL in Python

texts:
I want to convert ISIS cub into a tiff format. Following code the error
>> Array[Array == NDV] = np.nan
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: cannot convert float NaN to integer

#########################################

Code is here
import numpy as np
import gdal
from gdalconst import *
from osgeo import osr

# Function to read the original file's projection:
def GetGeoInfo(FileName):
    SourceDS = gdal.Open(FileName, GA_ReadOnly)
    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)
    return NDV, xsize, ysize, GeoT, Projection, DataType

# Function to write a new file.
def CreateGeoTiff(Name, Array, driver, NDV,
                  xsize, ysize, GeoT, Projection, DataType):
    if DataType == 'Float32':
        DataType = gdal.GDT_Float32
    NewFileName = Name+'.tif'
    # Set nans to the original No Data Value
    Array[np.isnan(Array)] = NDV
    # Set up the dataset
    DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType )
            # the '1' is for band 1.
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    # Write the array
    DataSet.GetRasterBand(1).WriteArray( Array )
    DataSet.GetRasterBand(1).SetNoDataValue(NDV)
    return NewFileName

# Open the original file
FileName = 'I29955002trim.cub'    # This is the ISIS3 cube file
                                  # It's an infra-red photograph
                                  # taken by the 2001 Mars Odyssey orbiter.
DataSet = gdal.Open(FileName, GA_ReadOnly)
# Get the first (and only) band.
Band = DataSet.GetRasterBand(1)
# Open as an array.
Array = Band.ReadAsArray()
# Get the No Data Value
NDV = Band.GetNoDataValue()
# Convert No Data Points to nans
Array[Array == NDV] = np.nan

# Now I do some processing on Array, it's pretty complex
# but for this example I'll just add 20 to each pixel.
NewArray = Array + 20  # If only it were that easy

# Now I'm ready to save the new file, in the meantime I have
# closed the original, so I reopen it to get the projection
# information...
NDV, xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName)

# Set up the GTiff driver
driver = gdal.GetDriverByName('GTiff')

# Now turn the array into a GTiff.
NewFileName = CreateGeoTiff('I29955002trim', NewArray, driver, NDV,
                            xsize, ysize, GeoT, Projection, DataType)


texts: GeoWebCache: Layer is not served by a WMS backend

texts:
I try to get legends for my layers in a little GIS application while using GeoWebCache.
Since I turned GeoWebCache in GeoServer on I do not get the legend image for my layers anymore.
I looked around and found this thread here: How to get WMS LegendGraphics using Geoserver and GeoWebCache?. However, my attempt to set the 'legendURL' by using '...LAYER" + map.layers[6].params["LAYERS"])' fails. The only thing I am getting back is that 'undefined' is not an object (evaluating 'mapPanel.layers[6].params').

Can somebody help me here? The application is temporarily here available.

Thanks a lot,
Thomas


texts:
I am using pandas, and folium to replicate the following example, but with my CSV data to produce a heat map.

Here's what my data looks like:

and here is the code I am using currently:

import pandas as pd
import numpy as np
import folium.plugins as plugins
from datetime import datetime, timedelta

    data = [(df3).tolist() for i in range(100)]

    m = folium.Map([40., -73.], tiles='stamentoner', zoom_start=6)

    hm = plugins.HeatMapWithTime(
        data,
        index=time_index,
        auto_play=True
    )

    hm.add_to(m)
    hm.save('map.html')

The example works, but I cannot input my own data, as it results in errors. Would anyone know how to plot the X, and Y data from my CSV or an approach to get me started to produce a heatmap?


texts: Getting Layername from selected feature using ArcObjects?

texts:
I have a IEnumFeature list of selected Features obtained through IApplication.Document.ActiveView.Selection.

From this selection I want to delete some feature, more specifically the features from one layer. Iteration through the IEnumFeature list is no problem but i cant figure out how to get the layername of each IFeature.

How can I obtain the layername from a IFeature?

I know how to get the ITable object from each IFeature but there seems to be no property that holds the name of the table, wich would suffice for my purpose. Because is know witch layer is created from wich table.


texts: Trying to count points in polygon and add a column to the table

texts:
I have a data variable called "schools" with the following data.

| School Name | Address | geometry              |
|---------------------------------------------- |
| School A    | 123 abc | c(123.321, 456,654)   |
| School B    | 321 cba | c(222.111, 333,444)   |
| School C    | 111 aaa | c(111.111, 222.222)   |
|-----------------------------------------------|

and the next variable is called "student_population"

| Region  | population | geometry |
|------------------------------------------|
| central | 2 | list(list(c(122.111, 333.111, 212, 321))) |
| West | 3 | list(list(c(122.111, 333.111, 212, 321, 112.331, 222.232))) |
| East | 2 | list(list(c(122.111, 333.111, 212, 321))) |

I want to count the number of schools in each region. hence, in the end i would get a table of the following

| Region  | population | number of schools in region | geometry |
|-----------------|-------------------------|
| central | 2 | 10 | list(list(c(122.111, 333.111, 212, 321))) |
| West | 3 | 20 |list(list(c(122.111, 333.111, 212, 321, 112.331, 222.232))) |
| East | 2 | 10 | list(list(c(122.111, 333.111, 212, 321))) |

I hope I am clear in explaining what I have done and what I want. I am having trouble solving this. In the end, I wish to tm_shape it. Tried using st_within but I do not know how to use it.


texts: Where is the forest field in openlayers plugin on Qgis?

texts: I am mapping out a region in my community and there is a vast forest area (I identified it from the openstreetmaponline map) and I will like to extract the shapefile polygon vector, but I don't know the 'key' value on the openlayers plugin on Qgis 2.14


texts:
In OpenLayers, I have a Vector layer with a StyleMap and SelectFeature and ModifyFeature controls, like this:

vector_layer = new OpenLayers.Layer.Vector('My Vector Layer', {extractAttributes: true,
            styleMap: new OpenLayers.StyleMap({'default':{
                strokeColor: "${customColor}",
                strokeOpacity: 1,
                strokeWidth: 2,
                fillColor: "${customColor}",
                fillOpacity: 0.5
            }}) });
map.addLayer(vector_layer);

selectControl = new OpenLayers.Control.SelectFeature(vector_layer,
                {onSelect: onFeatureSelect, onUnselect: onFeatureUnselect});
map.addControl(selectControl);
selectControl.activate();
modify_control = new OpenLayers.Control.ModifyFeature(vector_layer);
map.addControl(modify_control);

When I tick a checkbox on the page, the ModifyFeature control is activated:

modify_control.activate();
modify_control.mode |= OpenLayers.Control.ModifyFeature.ROTATE;
modify_control.mode |= OpenLayers.Control.ModifyFeature.DRAG;

The problem is this: when I remove the StyleMap, ModifyFeature shows the handles to drag and rotate the feature. With the StyleMap, they are not visible.
Any suggestions how to fix this?


texts: How to get routes in the database that intersects circle?

texts:
I have routes (linestrings) recorded in the database. I want to create a query which returns the routes that have some commonality with a circle whose center and radius are passed as parameter.

I tried this way:

SELECT *
FROM   routes
WHERE  ST_Touches(geom, (ST_Buffer(ST_GeomFromText('POINT(-12.988542214345332 -38.44563066959381)'), 100, 8)))

Is it correct?


texts: Summing Area of polygons within vector layer within zone using QGIS?

texts:
I have one vector layer with multiple polygons (purple). I also have a number of other vector layers which are catchment zones (I also have a layer where these catchments are merged).

Is there any way I can sum/calculate the areas of the purple polygons within each of the individual catchment zones?

I need this to work out a percentage of the area taken up by the purple polygons within each catchment.


texts: What geodetic reference system is QGIS using when providing latitude and longitude values for points using the field calculator?

texts:
When working in a projected map (project properties have a projection chosen), QGIS calculates longitude ($x) and latitude ($y) values for each point (or polygon centroid) as decimal degrees. (This is when using the field calculator to obtain lat/long). What geodetic reference system is QGIS using to provide lat/long values when using the field calculator?

I have made a few trials with a shapefile in WGS 84 and NAD 83 and the decimal.degrees values are the same in both shapefiles. However, they shouldn't be have the same values past the 5th or 6th decimal (based on using an online conversion tool from WSG 84 to NAD 83: https://tagis.dep.wv.gov/convert/)

EDIT: clarify the steps I made to test decimal.degrees values in a project using a projected CRS (NAD 83 MTM8)

Create new point shapefile
Create random cluster of points
Save point shapefile with WGS 84 as CRS (points_WGS84.shp)
Save point shapefile with NAD 83 as CRS (points_NAD83.shp)
Load points_WGS84.shp in QGIS
Open attribute table, open field calculator, create new field (long) with 8 decimal points, enter $x in the calculator box [points_WGS84.shp file]
Open attribute table, open field calculator, create new field (lat) with 8 decimal points, enter $y in the calculator box [points_WGS84.shp file]
Load points_NAD83.shp in QGIS
Open attribute table, open field calculator, create new field (long) with 8 decimal points, enter $x in the calculator box [points_NAD83.shp file]
Open attribute table, open field calculator, create new field (lat) with 8 decimal points, enter $y in the calculator box [points_NAD83.shp file]
Open attribute tables for points_WGS84.shp and for points_NAD83.shp
Compare values between both tables --> conclusion: they are identical to the last decimal


texts: Is ArcGIS for server 10.4 incompatible with the Silverlight API with regards to dynamic map layers?

texts:
We are currently using the Silverlight API version 3.0 and whenever we try to initialize an ArcGISDynamicMapServiceLayer from a 10.4 server we are getting InvalidCastException Unable to cast object of type 'System.String' to type 'System.Collections.ICollection'.

I know the latest (and final) Silverlight API is version 3.3, but upgrading breaks a lot of our existing code, and before spending any time trying to fix our code to upgrade to version 3.3 I was wanting to know if anyone else is having the same problems when running version 3.3

I've seen that other people have had the same problem https://geonet.esri.com/thread/173642, but I don't know what API version they were using.

If this is a known issue for the Silverlight API 3.3, then we would rather just stick with our current API version and continue putting our focus on moving away from Silverlight entirely.


texts: GIS software with Bayesian Network capabilities

texts: I am looking for any GIS software that has Bayesian Network capabilities.  Ideally I'd be able to define grid layers as discrete or continuous inputs in to the network and then be able to produce layers based on the predictions of given nodes.  I note there are several academic papers that seem to use a variety of different GIS and Bayesian tools but there doesn't seem to be any mention of an integrated solution.  I'm starting to think that a solution doesn't exist and I'm going to hack my own solution together.


texts: openlayers 3 and geolocation without HTTPS

texts:
I am creating an application where user can locate himself on map. the code is working fine in local development environment, but when i deploy it over server(IIS) the code won't work. here is code :

<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="chrome=1">
<meta name="viewport" content="initial-scale=1.0, user-scalable=no, width=device-width">

<title>location</title>

  <button id='buttonx' text='click'>Locate</button>

After a lot of research i found this link which summarize that the location feature will work over HTTPS only.

https://sites.google.com/a/chromium.org/dev/Home/chromium-security/deprecating-powerful-features-on-insecure-origins

I just want to ask do we have any process to bypass the feature.