texts: How to sort a SimpleFeatureCollection using GeoTools?

texts:
I have a featurecollection:

SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(dbFeatureType);
SimpleFeatureCollection collection = FeatureCollections.newCollection();
while (mifReader.hasNext()){
  in = (SimpleFeatureImpl) mifReader.next();
  SimpleFeature dbFeature = featureBuilder.buildFeature(null);
  for(AttributeDescriptor ad : adList){
    String name = ad.getLocalName();
    Object attrValue;
    attrValue = in.getAttribute(name);
    if(attrValue instanceof String){
      String string3 = new String((attrValue.toString()).getBytes("ISO-8859-1"), "cp1251");
      if(name.trim().equalsIgnoreCase(kadname.trim())){
        dbFeature.setAttribute("kadnum", string3);
      }
    }
    if (attrValue instanceof Geometry){
        idn++;
        com.vividsolutions.jts.geom.Geometry geom = (Geometry) in.getAttribute(name);
        dbFeature.setAttribute("id", idn);
        System.out.println("after insrt="+idn);
        dbFeature.setAttribute("deleted", deleted);
        dbFeature.setAttribute("the_geom", geom);
        dbFeature.setAttribute("status_id", status_id);

    }
    collection.add(dbFeature);
    }
}

Its all okay. But when I put this collection in Data Base i get a reversed items. So i want to sort or reverse items in this collection. But how to do this?
I try to use this:

        collection.sort( new SortBy()
{
    @Override
            public PropertyName getPropertyName() { return YourPropertyNameImpl("id") ; }
    @Override
            public SortOrder getSortOrder() { return SortOrder.ASCENDING; }
});

But I don't understand what is YourPropertyNameImpl.

UPDATE

Okey. I try to put features from my collection to another collection:

    Iterator iterator = null;
    int count;
    iterator = collection.iterator();
    while( iterator.hasNext()){
        SimpleFeature feature = (SimpleFeature) iterator.next();
        coll.add(feature);
        System.out.println("Coll= "+coll.features().next().getProperties("id"));
    }

But get same result(. In this case i need to read collection from last element to first. Its possible?


texts: Distortion when exporting from QGIS to DXF?

texts:
I'm rather new in the GIS world, and I'm having a weird problem that I have no idea how to solve.

I received a GIS map and I wanted to export one of the layers as DXF; I selected the correspondent CRS (WGS 84 / UTM zone 47N for Myanmar) and I generated a DXF, but when I opened the file it was somehow distorted.. I mean that the Y dimension was right, but the X dimension had "shrinked". In the attached image you can see that the new polygon (blue) matches the original boundary (orange) on the top (it also matches on the bottom), but not on the left and right sides.

Any idea of what can be the reason of this?

I'm using AutoCAD 2016 to open the DXF file.


texts: How to load the layer style from Postgresql with Python?

texts:
I would like to write a python macro that upon opening my QGIS project file will reload the default style for each layer from the postgres 'layer_style' table. I have had a look at the documentation and can only work out how this is done from a sqlite db or .qml file.

Any pointers would be appreciated.


texts: Ordnance Survey osgb:BNG spatial reference system not recognized by PostGIS

texts:
I'm trying to load an Ordnance Survey GML file into PostGIS using the ST_GeomFromGML function. The problem is that the spatial reference system that is used in the file, osgb:BNG, is not recognized by PostGIS:

# SELECT ST_GeomFromGML('<gml:Point srsName="osgb:BNG">
         <gml:coordinates>527347.0,181346.0</gml:coordinates></gml:Point>');
ERROR:  unknown spatial reference system

Is there any workaround? Does "osgb:BNG" have an alternative name that is recognized by PostGIS? Thank you.

EDIT: Thanks to everybody for their answers. I ended up replacing "osgb:BNG" with "EPSG:27700" in the GML file (in the srsName parameter), although the other two options (modifying the spatial_ref_sys table, or using the ST_GeomFromGML function http://postgis.net/docs/ST_GeomFromGML.html with two arguments for PostGIS 2.0+) are perfectly valid too.


texts: Split line at points where other lines touch it in PostGIS

texts:
I am trying to create an application that will require a network to be rebuilt several times in an iterative process, adding connections as it progresses.  Every iteration needs to use the network created in the last iteration to create the next network.

I am using pgrouting to calculate routes along the network during every iteration, and hoping that building a pgrouting topology does not try to node every intersection, because there are several under/overpasses in the data that I would prefer to simply leave as they are.

However I need to node the new connections that are being made to the existing network, so that I can include them when rebuilding the next network iteration.  I have been able to draw a line connecting new points to the network, as well as creating a set of multipoints where the line meets the existing network.  I want to split the existing network only at these points, but have all segments retain the attributes of the original line.

Essentially, this is the situation:

I want to manually split the pink line at all the pink dots, so that I can rebuild a pgrouting topology, and I want all the subsequent segments to retain the attributes of the pink line itself.

This is the code that I am using.  I just keep receiving the original line's geometry, un-split.

SELECT ST_AsEWKT(ST_Split(N, P)) AS geom
FROM (SELECT
  starting_network.geom as N,
  cut_points_multi.geom as P
  FROM starting_network, cut_points_multi) AS foo;


texts: ERROR 99999 during arcpy.Project_management

texts:
I have a script that I run in ArcMap 10.5.1 that iterates through feature classes in a file geodatabase and projects each one to a WGS84 shapefile.  This works smoothly for most feature classes in the geodatabase, but there is one feature class that is giving an odd error:

Traceback (most recent call last):
File "E:\ProvingGrounds\SCAT\Script\SCAT Shapefile Export 02.py", line 38, in <module>
    arcpy.Project_management(fc, outShape, spatialRef)
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 9523, in Project
raise e
ExecuteError: ERROR 999999: Error executing function.
The input featureclass fields and the Output fields do not have the same number of fields
Failed to execute (Project).

This seems strange since arcpy.Project_management is creating the output from the input.  The feature class in question is as simple as it gets, a single feature polyline with no attributes other than ObjectID and Shape_length.

Here is the script that generates the error:

import arcpy, os
arcpy.env.overwrite = True

scatFolder = arcpy.GetParameterAsText(0)
gdb = os.path.join(scatFolder, "scatprocessing.gdb")
shapeFolder = os.path.join(scatFolder, "ShapeFiles")
spatialRef = arcpy.SpatialReference("WGS 1984")

if arcpy.Exists(gdb):

    if not os.path.exists(shapeFolder):
        arcpy.AddMessage("\nCreating shapefile folder\n")
        os.makedirs(shapeFolder)

    arcpy.env.workspace = gdb
    featureClasses = arcpy.ListFeatureClasses()
    for fc in featureClasses:
        arcpy.AddMessage("Exporting " + fc + " to shapefile")
        outShape = os.path.join(shapeFolder, fc)

        arcpy.Project_management(fc, outShape, spatialRef)

The standalone Project tool in the Management Toolbox works as advertised on the same feature class using the same parameters. It is only when iterating through feature classes (to do the same thing) that it throws the error on that feature class.


texts: Error clipping raster using a mask layer

texts:
I am trying to clip a raster with a polygon with QGIS 3.0 on High Sierra OSX (I've done it before with OSX Sierra and QGIS 2.8, no issues), however I keep getting an error (see inputs in screenshot below)

Here's the error:

*Processing algorithm…
AlgorithmDialog
Input parameters:
{ 'INPUT' : '/Users/felipearrospide/Dropbox/WATS/W04-18 Natura/2- Trabajo/gt30w100s10_dem/gt30w100s10.dem', 'MASK' : '/Users/felipearrospide/MEGA/PROYECTOS/01 2018/W04-18 Natura/2- Trabajo/Cuenca PFV Itahue.shp', 'NODATA' : 0, 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'KEEP_RESOLUTION' : False, 'OPTIONS' : '', 'DATA_TYPE' : 5, 'OUTPUT' : '/Users/felipearrospide/MEGA/PROYECTOS/01 2018/W04-18 Natura/2- Trabajo/DEM/trial.tif' }
GDAL command:
gdalwarp -ot Float32 -of GTiff -cutline "/Users/felipearrospide/MEGA/PROYECTOS/01 2018/W04-18 Natura/2- Trabajo/Cuenca PFV Itahue.shp" -crop_to_cutline "/Users/felipearrospide/Dropbox/WATS/W04-18 Natura/2- Trabajo/gt30w100s10_dem/gt30w100s10.dem" "/Users/felipearrospide/MEGA/PROYECTOS/01 2018/W04-18 Natura/2- Trabajo/DEM/trial.tif"
GDAL command output:
/bin/sh: gdalwarp: command not found
Execution completed in 0.09 seconds
Results:
{'OUTPUT': <QgsProcessingOutputLayerDefinition {'sink':/Users/felipearrospide/MEGA/PROYECTOS/01 2018/W04-18 Natura/2- Trabajo/DEM/trial.tif, 'createOptions': {'fileEncoding': 'System'}}>}*

Am I missing something? Should I install/reinstall/uninstall everything?
I am truly lost :(


texts: How to add TomTom's traffic flow map images to Leaflet?

texts:
I have a raster image from TomTom which shows traffic flow. The image is created after sending lat, long coordinates, that is one point to the GET serves with the zoom value, then I get the image title.

How to lay this image on the map in the right position?

Note: There is no Bounds available.

If I use L.tileLayer, I get repeated image see image:


texts: rasterize gradient from shapely object with rasterio

texts:
I have an array of horizontal and non-horizontal 3d-polygons as shapely objects, and I want to rasterize them into a depth map

For horizontal ones, it's just drawing their height as color, but I also need to render a gradient for the non-horizontal ones; how should I approach that?


texts:
I'm trying to look at some landsat images and thought I would give the RS&GIS a go, as a way of creating true and false colour composites. I'd like to visually inspect the vegetation my project area.

Unfortunately, when I try to run the plugin it produces no outputs. The dialogue box notes the following "IndexError: list index out of range".

I've already tried the below but without success:

removing all spaces from the file path
using different raw inputs (compressed and uncompressed)
with and without the clipping function


texts: Using the excludelayer tag in Flex 2.3.1 with a cached service

texts:
I'm working on our first Flex server page using the Flex Viewer 2.3.1. The Legend widget works, but I cannot get the excludelayer tag to work. Here's where the problem probably is - it's a cached map service. Can the excludelayer tag be used on a cached map service? The layer I'm trying to exclude is an image service. Thanks in advance.

Here's the snippet of code:

<?xml version="1.0" ?>
<configuration>
    <usebasemaps>true</usebasemaps>
    <useoperationallayers>false</useoperationallayers>
    <respectcurrentmapscale>true</respectcurrentmapscale>
    <excludelayer>SJCounty_2009_1ft</excludelayer>
</configuration>

The image service is within the base map layers. I'm trying to remove the entire image service from the legend, and not just the boundary/footprint,etc.


texts: Installing ArcGIS Desktop on Windows 7 (64 bit)

texts:
I just reformatted my PC with Win7 64bit to try and take advantage of 64-bit python processing.   However I'm not sure what to do here:
Our office runs ArcGIS 10.1, TOAD and SDE for Oracle.     Toad says I need to install a 64-bit client before I can connect to Oracle.   The ArcGIS Resource Center says before I can connect to SDE I need to install a 32-bit oracle client.    I am confused.


texts: PostGIS query to retrieve the largest polygon for multi-polygons by grouping on attribute field

texts:
I'm working with the Natural Earth 10m Admin 1 States Provinces polygon data. This dataset contains multi-polygons for each state / province in the entire world. What I'd like to do is to grab the largest polygon for each state / province, then create a centroid that is inside that polygon for labeling purposes in web cartography.

I went ahead and calculated the area for each polygon using ST_Area(geom::geography) for the dataset and so far my PostGIS query looks like this, which returns the the largest polygon for a single state /province:

SELECT area, name, ST_Centroid(geom) geom
FROM ne_10m_admin_1_states_provinces
WHERE area IN (
      SELECT area
      FROM ne_10m_admin_1_states_provinces
      WHERE name = 'state province name'
      ORDER BY area DESC
      LIMIT 1
)

However I'm unsure of how to determine the best way to iterate over all features in the dataset to grab the polygon for each state / province with the largest area. Would the best way to do this be to create a function in PGSQL? Or is there a simpler way via plain old SQL?


texts: ArcSDE on SQL Server truncating fields even after changing Field Length

texts:
I created a SQL Server Enterprise Geodatabase for our Fiber Network. Sometimes our strand descriptions get pretty long. I ran into an issue where some of the descriptions were getting cut off at 254 which is what I had set the length of the description field originally. I made sure everyone was out of the database and I changed the field length to a sufficent 1000 characters in SQL Management Studio.

The issue is that it will accept the additional characters UNTIL the user saves the edits, then it will truncate back down to 254 characters. I have been testing and testing and I compressed the database and I cannot figure out what the issue is.


texts:
How can I resolve problematic kerning in GeoServer labels?

Example:

Problematic kerning is in "San Glissando"

Another example:

Geoserver documentation contains only:

<Font>
       <CssParameter name="font-family">Arial</CssParameter>
       <CssParameter name="font-size">12</CssParameter>
       <CssParameter name="font-style">normal</CssParameter>
       <CssParameter name="font-weight">bold</CssParameter>
</Font>

I would like to find something like the CSS parameters, which can change spacing between letters in labels.


texts:
Is there an implementation of Openstreetmap that closely mimics the appearance of 1:25k OS Explorer or 1:50k OS Landranger maps? By "mimicking the appearnace" I mean same colour scheme, same POI symbols, same typography and labeling conventions, contour lines, etc.
If the answer to 1. a is "no", is it possible to create such an implementation? What might be a good resource to learn how to do so?

Like many Brits, I grew up with paper OS maps and am still more comfortable reading them than a style such as that used at http://www.opentopomap.org/


texts: If condition for a table in geodatabase in ArcPy?

texts:
I want to do what would be in a dataframe:

if df['col'] = 1:
    do this..

but for a table in a geodatabase:

if col = 1  in arcpy.ListTables("qt"):
   do this..

Purpose:

For the rows that have in column 'col' the value 1 in a table in a geodatabase  do a piece of code.

How to do this in ArcPy?

Update

sample table:

names  col
A       0
B       1
C       1
D       0
E       1

the block of code after the condition check, therefore should be implemented in rows B,C,E.


texts: How to check if a rotated rectangle is contained in another rotated rectangle?

texts:
I need to check if a rotated rectangle A is inside a rotated rectangle B. That is, if B is completely contained in A. I have the following information stored in a rotated rectangle: the four coordinates and its angle.

Is there is an efficient method to do this? I know that the Separating Axis Theorem can check a intersection situation, but is this method still valid to check a containment situation?


texts:
I downloaded Landsat 8 NDVI bands from the USGS (instead of calculating them myself). The QA band that came with has values different than I see in the Landsat 8 user guide.

Using R I did a quick sort(unique(as.vector(image[[1]]))) and got these values:

322 324 328 336 352 368 386 388 392 400 416 432 480 834 840 848 864 880 898 904 912 928 944 992

Here is the actual QA raster. Guessing high values are clouds.


texts: How to disable layer(not hide) in OpenLayers 3

texts:
In my map, I have  multiple layers.
For certain operations I want to only work on one layer. The other layers should be visible, but disabled, so that they can not be selected or interacted with in any way.


texts: Populate empty attribute table from constantly growing CSV in QGIS

texts:
I have a Python script running that will generate a line feature from a point depending on a degree of rotation around that point. I have a non-spatial CSV setup to get fed this rotation data from another source. The CSV will be constantly updating and generating new data, and will generate new lines in QGIS as they arrive.

As you might have inferred, I want this CSV to effectively act as a dynamic attribute table for this line feature class, but I'm having trouble figuring out the right way to perform a table join on the line feature. I have an "ID" and a "LOB" (rotation) field set up in both the CSV and the empty attribute table, so my expectation was that the tables would join and the CSV data would carry over and appear in the table. But this is not the case, as I somehow need to get the number of rows in both "ID" columns to match prior to the join. But this doesn't seem like a solution, as the number of columns will be rapidly changing.

Are there any reasonable solutions to this dilemma? This is very easy to do with spatial text-delimited point layers, but mind boggling with non-spatial data for lines.

Below is a screenshot for reference. The top table is the CSV, the bottom is the attribute table it's "joined" to. I need the number of rows in the attribute table to match that of the CSV.


texts: Creating Dashboard Templates in ArcGIS Online

texts: I have a dashboard that took some time to build and I would like to use it as a template for other projects. Is this possible?


texts: As a land surveyor, I often bring cad data into a GIS and create polygon representing parcels.  Traditional tools I have for cleaning topology are inappropriate so I am looking for something different - even if it means teaming up with someone to build additional extensions.  (I scanned the list of existing extensions and some pieces are already there).  My idea is to fix minor topology errors by constraining key geometric elements. That is, straight polylines and circular arcs which share a common centre point.  A key requirement is support for true circular curves.  I could keep the data in SQL Server but prefer a GUI for a better user experience.  As I am not familiar with QGIS, I am looking for some indication if the data structure used by QGIS can support this.


texts: How to apply texture on a 3D model using drone images with python?

texts: I have generated a 3D model in .OBJ format using UAV (drone) images. But this model doesnt have texture so i want to have it textured using the original UAV images. And this should be done in python so that I can include in my original code.


texts: removing spurious lines in raster before or after vector conversion?

texts:
I'm looking for a feature/function in either python opencv/ndimage/gdal/etc. or qgis/saga/grass that will help clean up a skeletonized image that I'd like to vectorize into a shapefile of lines/polylines.

I spent considerable effort producing this network and I'd like to vectorize it and remove 'jaggies' (short lines to go nowhere) and also take the shortest path through the clumps of webby areas to produce as straight lines as possible.

I realize this is a multi-step process and there may be possibilities in the python space or in gis.  I know the steps I would do if I did it from scratch but it would take a while!


texts: Convert PNG raster MBTIles tiles to JPG raster MBTIles Tiles

texts:
I have a set of raster png MBTiles that are way too big in terms of file size. I need to get them down to close to half. I was hoping to achieve it by converting them to JPG Mbtiles, as they are currently PNG.

Do you know any tools that can achieve that?


texts: How to calculate a field based on other attributes

texts:
I am trying to update a field "SZone" with the information from field "Zone", based on field "PointOrder". I am trying to use Calculate Field and the following code, but it is not updating.

I am using the expression:

new(!Update!)

def new(Update):
  if (PointOrder ==0):
    return Zone
  else:
    return ""


texts: I'm using eCognition to mapping roof materials from SPOT 6 image. I have used nearest neighbour classification but the result is not so good. I'm trying to apply a random forest classifier if it is available in eCognition (maybe need to plug in).


texts:
I would like to analyze NDVI in the region of Bangkok.

Here is the code :

{
    //Define time range
    var startyear = 2006;
    var endyear = 2016;

    //Set date in ee date format
    var startdate = ee.Date.fromYMD(startyear, 1, 1);
    var enddate = ee.Date.fromYMD(endyear, 12, 31);

    ///Import a fusion table of your study area
    // The Ca basin
    var Ca = ee.FeatureCollection(geometry)

    //filter on date and area (bounds)
    var l5images = l5.filterDate(startdate,enddate).filterBounds( Ca);
    var l7images = l7.filterDate(startdate,enddate).filterBounds( Ca);
    var l8images = l8.filterDate(startdate,enddate).filterBounds( Ca);

    ///Include a funtion to remove the clouds
    //set cloud threshold
    var cloud_threshold = 40;

    // Select the red and NIR bands
    l5images = l5images.select(["B4","B3"]);
    l7images = l7images.select(["B4","B3"]);
    l8images = l8images.select(["B5","B4"]);

    // calculate ndvi from landsat 8
    function l8ndvi(img) {
      var ndvi = img.normalizedDifference(['B5', 'B4']).rename('NDVI');
      return img.addBands(ndvi);
    }

    // calculate ndvi from landsat 5 & 7
    function l57ndvi(img) {
      var ndvi = img.normalizedDifference(['B4', 'B3']).rename('NDVI');
      return img.addBands(ndvi);
    }

    // calculate ndwi for each image in imagecollection
    var l5ndvi = l5images.map(l57ndvi);
    var l7ndvi = l7images.map(l57ndvi);
    var l8ndvi = l8images.map(l8ndvi);

    // combine all data in single image collection
    var allcollection =  ee.ImageCollection((l5ndvi.merge(l7ndvi)).merge(l8ndvi));

    // add map to canvas
    Map.addLayer(allcollection);
    Map.centerObject(allcollection,8);
    print()
}

Google Earth Engine Link that I used

But then there is trouble, but I don't really understand what is it about.

So Is there anything I can do about this. I am still very new to programming.


texts: Programatically create WMS and specify data sources in QGIS Server from a web app

texts:
Is it possible to programatically create a WMS and specify data sources in QGIS Server from a custom web app?

Let's say I have a QGIS project being served up by QGIS Server and accessible via a web client (Lizmap or QWC2 for example). Features are stored in PostGIS (or MSSQL or GPKG). Symbols and labels are all set up and looking good:

http://localhost/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&MAP=/home/ubuntu/qwc2/proj1.qgs

Now suppose I want to create a second WMS that points to a new set of similar PostGIS layers located nearby.  Rather than manually create a second QGIS project with layers renamed and data sources changed:

http://localhost/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&MAP=/home/ubuntu/qwc2/proj2.qgs

Is there a way to programatically create this new QGIS project (new WMS) from say a custom web app using Python?  The only difference between proj1.qgs and proj2.qgs would be layer names, data sources and the bounding box.  The underlying PostGIS layer schema (geometry, fields) are identical.

I have a QWC2 web map that contains multiple themes (one QGIS project per theme).

My ultimate goal is to create a new theme by specifying a new data source in a custom web app.  The new WMS would then get created automatically.


texts: Adjusting QGIS Print Composer output resolution?

texts:
I am using QGIS 2.12.

I would like to adjust the Print Composer's output resolution for "Export as Image" and "Export as PDF".

The default is 300 dpi, but I would like 600 or 1200.

Multiple websites say that it can be done, but they don't say how.


texts: Add OpenLayers3 controls to qgis2web project

texts: I am trying to get some of the control sets offered by OpenLayers into a project exported using QGIS2WEB. I am looking to use the export to Pdf function and have successfully created a test page using it. However I do not know how to integrate the function/script into the qgis2web project. Eventually I would also like to also add the draw feature control


texts: How to measure distance between points based on elevation?

texts:
I need to be able to measure distances between points, however the distance need to be calculated in relations to elevation. The points are homes of storytellers from the 19th century and the places which are mentioned in their stories. The distance therefore has to be "walking distance". A path along a valley will likely be shorter than a path over a mountain even though the actual straight distance is shorter. Attached is a screenshot illustrating my thinking. In the image, paths A and C would therefore be calculated shorter than path B.

The points are from a CSV file but I also have a raster layer with the elevation data.


texts: How to replace data in parts of raster A with the equivalent data in raster B?

texts:
I have the following problem: I am currently working with Landsat 5 imagery from the tropics including quiet high cloud cover. I have classified raster images with the different land use categories I need and "Clouds" as well as "Cloud shadow" as categories where there appear on the image.

I now extracted the "Clouds" and "Cloud shadow" parts of my classified image (Extraction toolset --> Extract by Attributes) and converted these to polygons. Now I have a layer with polygons for every cloud and every shadow of my classified image. I was hoping that I can now take a different classified image from the same year with less clouds in the areas where my original image has a lot of clouds and extract the data from this new classified image that falls within the polygons that I created. Then I wanted to take this new data and "overlay" my original classified image with the extracted data that falls within the area that is covered by clouds. The result would be that I replaced the areas covered by clouds with the data from the other classified image of the exact same area.

However, I can not find any tools to do this overlay-like operation.
If someone knows how to extract data from a classified Raster image using polygon files (similar to a mask) and how to replace the data from the original classified raster image with the newly obtained data I would highly appreciate your help.


texts: Png, prj, pgw same name, but png doesn't show up after loading. What SRC to take?

texts:
someone gave me a dataset with the three types of files, but when I load as a raster the png.files it doesn't show up on the maps. They ask me what SRC to take and I don't really know which one to choose, I tried a bunch (the one of the others layers, one that seems on my part of the world etc..) none seems to work. As I have the proper pjr and pgw files that go with it, why it doesn't read it as the geographic reference. .. Tried several things and read some to, but I still don't understand.

This is the content of the PGW files:

1.0000
0.0000
0.0000
-1.0000
240591.5000
4991759.5000

Prj file:

PROJCS["NAD83_/_MTM_zone_8",GEOGCS["NAD83_/_Géographique",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.2572221010,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0.0],UNIT["Decimal_Degree",0.01745329251994330]],PROJECTION["Transverse_Mercator"],PARAMETER["Latitude_Of_Center",0.0],PARAMETER["Longitude_Of_Origin",-73.5000000000],PARAMETER["Scale_Factor",0.9999000000],PARAMETER["False_Easting",304800.0],PARAMETER["False_Northing",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG","32188"]]


texts: Converting a NAD27 Lambert Conformal Conic custom projection to WGS84 using QGIS

texts:
I am trying to convert some points from NAD27 Lambert Conformal Conic (custom projection) to WGS84 using QGIS 2.8.9. I am using the datum transformation parameters for El Salvador published here: http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

I have created a custom projection using the following proj4 string:

+proj=lcc +lat_1=13.31666666 +lat_2=14.25 +lat_0=13.78333333333333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +k_0=0.999967040 +datum=NAD27 +ellps=clrk66 +units=m +no_defs +towgs84=0,125,194,0,0,0,0

It seems the +towgs84 parameters are not being applied. If I drop the +datum=NAD27 it seems it works. I have seen this suggestion in On the fly transformation From WGS84 to NAD27 with proj4 code, but I don't understand why it works. According to this in Reprojecting between NAD27 and WGS 84 part III, the +datum=nad27 parameter overrides the +towgs84 parameter, so I can use only one of them.

So my question is, can I specify a datum and some transformation parameters to WGS84 in proj4, or I just need to specify the ellipsoid? (Although it's an ellipsoid, not a datum)


texts: Intersect tool of ArcGis gives multiple intersecting points

texts:
I am using .net in arcGis 10.0 to develop tools.I have two polyline feature classes. I want to get intersected point of this 2 feature classes. for this purpose, I have used INTERSECT tool of arcmap. I want Just one intersected point of 2 features. But INTERSECT tool gives 4 points. How can I get one intersected point of two line feature ? Can I use another tool for this purpose??

I have an attached image .In the image, selected 4 feature denotes points in the intersected points. But I want one points.


texts: Exporting to 2D (Arcview 3.0 compatible) shapefile in QGIS?

texts:
I need to run my shapefiles through a web mapping application which is incompatible with 3D shapefiles.

At present I run them through the FME Universal Translator in Mapinfo, which enables you to specify shapefile type.

I'm trying to find equivalent functionality in QGIS but haven't been able to so far.

Anyone know if this is possible?


texts: OpenLayers 4 get style from GeoJSON not find

texts:
We created an API that returns a GeoJSON in accordance with this styling guide but then OpenLayers doesn't find the style key once I want to style each feature, here my code

var appLayer = new ol.layer.Vector({
  source: appSource,
  style: function (feature, resolution) {
    console.log(feature.getProperties()) // <== all geojson properties
     return [new ol.style.Style({
       fill: new ol.style.Fill({
         color: feature.get('...')
       }),
       stroke: new ol.style.Stroke({
         color: feature.get('...'),
         width: feature.get('...')
       })
    })]
  }
})

and here my GeoJSON file

What I'm doing is creating a Vector layer and then add the features from the GeoJSON as features.

$.ajax({
  type: 'GET',
  url: 'https://www.traffwebdev.uk/tro/v1/restrictions?epsg=&x=383666+&y=398197&area=&radius=50&output=',
  dataType: 'json',
  success: function (data) {
    // reads and converts GeoJSon to Feature Object
    var features = (new ol.format.GeoJSON()).readFeatures(data, {
      dataProjection: 'EPSG:4326',
      featureProjection: 'EPSG:27700'
    })
    appSource.addFeatures(features)
  }
}).then(() => appLayer.getSource().refresh())


texts:
I have gridded raster files with 6 km resolution, and each cell has its own respective Snow Cover value. Now, I want to change the resolution of this raster to 500 m, and based on the 6 km pixel values assign new Snow Cover values to each pixel of 500 m resolution.
Which method in ArcGIS Desktop 10.7 (interpolation or resampling technique) would best suit this kind of conversion and data?


texts: I'd like to know what the difference is between place_id and osm_id. I understand osm_id is some sort of reference number + more convenient for OSM to find location v. lat/lng. If this is the case, what is place_id and why, in nominatim.openstreetmap.org, do you query by that param instead? Perhaps it's internal to Nominatim? Can someone explain the differences, the only result Google spat out was in German and a bit of a mess when translated (http://forum.openstreetmap.org/viewtopic.php?id=19671).


texts: I need to merge two different layers but they don't have the same geometry, in particular one contains points and one is a polygon (province). Is it possible to do so by converting them to the same geometry and what is the best way to do it?


texts: How disable login prompt when asking for featuretypes?

texts: Is it possible somehow disable login prompt when requesting  http://xxxxxxxx/geoserver/rest/workspaces/WORKSPACE_NAME/featuretypes.json?


texts: Updating arcpy.Array() syntax from Python 2.7 to 3.6.6 for ArcGIS Pro 2.3?

texts:
I was working through an exercise designed to teach basic Python skills in ArcMap using Python 2.7. The file below consists of 34 pairs of x,y coordinates. The script creates a new empty polyline feature class, and then use this list of coordinates to create a new polyline. This requires the use of an insert cursor to create new rows and an array object to contain the point objects. Here is the original proposed script:

In the Spyder IDE with the ArcGIS Pro conda installation for 2.3, I hit a number of exceptions that I worked through as follows:

import arcpy
import fileinput
import string
import os
from arcpy import env
env.workspace = "C:/Workspace"
env.overwriteOutput = True
arcpy.CreateFeatureclass_management(r"C:/Results", "newfc", "POLYLINE", None, "DISABLED", "DISABLED", None, None, 0, 0, 0, None)
infile = "C:/coordinates.txt"
cursor = arcpy.da.InsertCursor("newfc.shp", ["SHAPE@"])
array = arcpy.Array()
for line in fileinput.FileInput(infile):
    ID, X, Y = (float(infile.split(line," "))
    array.add(arcpy.Point(X, Y))
cursor.insertRow([arcpy.Polyline(array)])
fileinput.close()
del cursor

However, there is a syntax error related to the arcpy.Array() method.

      File "C:/create.py", line 21
array.add(arcpy.Point(X, Y))
    ^
SyntaxError: invalid syntax

Here is the point file.

My knowledge is limited and wondered if someone knows if this is a Python version issue and where to find more info, or if they have an answer to the syntax error problem I articulated above.


texts: How to get the intersected region(s) of multiple polygon layers

texts:
I have multiple layers with polygons.

I want the regions that are intersected by all the polygons. The shared, common regions.

An example use case would be finding the shared area of multiple raster files if you had their footprints as polygons.

How can I achieve that with built-in QGIS tools?


texts: Assign symbology to multiple layers in one operation?

texts:
Is it possible to assign the same symbol to multiple layers of the same geometry type within an ArcMap document in one operation?

For example, I was hoping to:

Select multiple layers in the TOC
Click on the symbol "patch" or right-click an get to the "Symbol Selector".
Apply the same symbol to all selected layers.

Obviously that doesn't work, or I wouldn't be asking.  I found an old 2003 posting that talks about how to do this with VBA code.  Can this type of operation be done from the GUI in 10.x, or am I just out of luck?


texts: Calculating shape area for polygons in different files using Python

texts:
I need to calculate the shape area (and possible the perimiter) of polygons for over 300 files. I tried to use a Python code but I am not very familiar with Python so when it returns an error I do not know where the problem is. My Python code looks like this:

import arcpy
from arcpy import env
import os
import csv

outputFilePath = 'outputPath'
sourceFolderPath = 'inputPath'

for root, dirs, files in arcpy.da.Walk(sourceFolderPath):
    for name in files:
        if name[-4:] == '.shp':
            arcpy.env.workspace = os.path.join(root, name)
            arcpy.AddField_management(name[:-4], "Shape_area", "DOUBLE")
            arcpy.CalculateField_management(name[:-4], outputFilePath + '/' + name[:-4], "Shape_area", "PYTHON_9.3")

I received this error message:

Traceback (most recent call last):
  File "<string>", line 15, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 3290, in CalculateField
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 3287, in CalculateField
    retval = convertArcObjectToPythonObject(gp.CalculateField_management(*gp_fixargs((in_table, field, expression, expression_type, code_block), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.

Any suggestions?


texts: Summing of individual image's statistic over ImageCollection Google Earth Engine

texts:
I have an ImageCollection containing FIRMS fire image in 2014. I want to count the number of fires in the year for a predefined region. I can find a way to count the number of fires for an individual image through the use of ee.Image.reduceRegion.
However, I cannot figure out a way to summing up individual image's fire count to derive total fire count for 2014.
Here's my code:
        // Load Vietnam shapefile
    var vn_shape = ee.FeatureCollection('users/chucmd/VNM_adm0');
    print (vn_shape);

    // Load fire counts image
    var fire = ee.ImageCollection('FIRMS').filterBounds(vn_shape).filterDate('2014-01-01', '2014-03-30');

    // Clip image to Vietnam
    var fire_vn = fire.map(function(image) {return image.clip(vn_shape);});
    print (fire_vn);

    // Filter fire with more than 50% confidence and add a new band representing areas where confidence of fire > 50%
    var filterConfidence = function(image) {
      var line_number = image.select('line_number');
      var confidence = image.select('confidence');
      var conf_50 = confidence.gt(50).rename('confidence_50');
      var count_band = line_number.multiply(conf_50).rename('count');
      return image.addBands(count_band);
    };
    var fire_conf_vn = fire_vn.map(filterConfidence);
    print (fire_conf_vn);

    // fire count for individual image
    var img1 = ee.Image(fire_conf_vn.first());
    print (img1);
    var countObject = img1.reduceRegion({
        reducer: ee.Reducer.countDistinct(),
        scale: 1000,
        geometry: vn_shape
      });
    print (countObject);

And the result:


texts: How to burn coordinate point of polyline to the pixel value while converting it from vector to raster

texts:
I am trying to plot a proximity raster for secondary roads in a particular area. What I have is a set of polyline representing the roads. Now, I want to create a proximity raster based on the polyline.
First of all I need to create a raster from the polyline which I am perfectly able to. Now, how can I insert burn in value of the pixel to be it's lat and long coordinate?


texts: Couldn't destroy QGIS Custom Edit Form Factories

texts:
I created custom form factories using this tutorial. However, the factory destructor QgsEditorWidgetRegistry(). However, it doesn't destroy them. I just used the above destructor at the beginning of the creation of factories like this but there is no effect.

QgsEditorWidgetRegistry() # Destroy existing factories
myFactory = WidgetWrapperFactory()
QgsEditorWidgetRegistry.instance().registerWidget( "MyRuleWidget", myFactory )

Is there anything that I am missing or is it a QGIS bug?


texts:
I am looking to change versions of an .sde through Arcpy to then afterwards,  run a few geoprocesses on it. I do not want it to automatically update after that.

So far I have managed to create a version through ArcPy, but I can't for the life of me figure out how to connect to it. The "Changeversion_management" only does this on 1 feature class and keeps throwing me errors such as the following:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Feature Layer or Table View: Dataset ZN_KGD@pkgd.sde does not exist or is not supported
Failed to execute (ChangeVersion).

And

ERROR 000800: The value is not a member of ZN_KGD.ExPimp_p9 | ZN_KGD.Tool_test
Failed to execute (ChangeVersion).

I would like to change the version of all the feature classes in the TOC, or at least connected to the same .sde. The script I am using to create and change the version is the following:

import arcpy

Project_number = arcpy.GetParameterAsText(0)

inWorkspace = "C:/Users/Lansr/AppData/Roaming/ESRI/Desktop10.3/ArcCatalog/ZN_KGD@pkgd.sde"
parentVersion = "ZN_KGD_DIENSTVERSIE"

# Execute CreateVersion
arcpy.CreateVersion_management(inWorkspace, parentVersion, "ExPimp_{}".format(Project_number), "PUBLIC")

arcpy.ChangeVersion_management('ZN_KGD.ZN_KGD_DIENSTVERSIE',
                               'TRANSACTIONAL',
                               "ExPimp_{}".format(Project_number),
                               '')

The feature layer and sde are not defined in the script, but are purely in the TOC.

It is important to note that the top error happens when I fill in the name of the .sde and the bottom error happens when I fill in the name of a feature layer. My knowledge of ArcPy is limited at best. I am using ArcGIS 10.3.1.


texts: Line intersection with circle on a sphere (globe or earth)

texts:
How can I find the intersection point of a line and a circle on a sphere (earth surface)
Please, if any one knows a good solution help.


texts:
I have a point layer, with angle(real), width(real) and height(real).

I want make a ellipse marker that take this values.
But how I can offset the marker to left border (see http://1drv.ms/1lErGUp) . I try but «anchor point» (Hleft or Hright) don't make it (Vup and Vdown is alright).

I try to make an expression in «data defined properties» - , something like :

  "width"/2 || ',0'

or

  "width"/2 || ',0.0'

etc.


texts:
I'm trying to create a pop-up based on a conditional IIF statement and am struggling with the semantics of arcade. This is my code. (Note: I've tried 10 both in and not in quotes)

IIf(IsEmpty($feature["PERCENT_A"]), "10", $feature["PERCENT_A"])

Once I base my pop-up on this expression, it looks like the false condition is working but the true condition is not; or else the popup is still just displaying my PERCENT_A field and applying no conditions to it.

I'm working in arcgis online and in the pop-up contents settings I've selected Display: 'A description from one field' and then selected my saved expression from the dropdown.


texts:
I am trying to to create multiple input autocomplete search fields that function same as the demo at https://www.mapbox.com/mapbox-gl-js/example/mapbox-gl-geocoder/  . However, I do not want the Search input box to be inside the Map box, I want to have the Search field outside of the map. I have not found a good example of how to create this using the mapbox gl examples and api docs.

I am looking for something similar to Google Autocomplete feature (https://developers.google.com/maps/documentation/javascript/examples/places-autocomplete) but in Mapbox GL where I can just reference and id (ie, id="from_address") and Mapbox will attach the autocomplete behavior.


texts:
I am trying to add a GeoTiff to GeoServer but it is appearing as Unknown in Native SRS. Any idea what could be causing this?

Details for projection are.

PROJCS["NAD_1983_StatePlane_California_VI_FIPS_0406_Feet",
  GEOGCS["NAD83",
  DATUM["North American Datum 1983",
      SPHEROID["GRS 1980", 6378137.0, 298.257222101, AUTHORITY["EPSG","7019"]],
      TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
      AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4269"]],
  PROJECTION["Lambert_Conformal_Conic_2SP"],
  PARAMETER["central_meridian", -116.25],
  PARAMETER["latitude_of_origin", 32.166666666666664],
  PARAMETER["standard_parallel_1", 33.88333333333333],
  PARAMETER["false_easting", 21527734.722222216],
  PARAMETER["false_northing", 5381933.680555556],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["standard_parallel_2", 32.78333333333333],
  UNIT["foot_survey_us", 0.30480060960121924],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH]]


texts: Two different attribute symbology in QGIS

texts: I have protected species map in which I have different attributes two of them are taxon group and common names. I want to display the symbol that within  each taxon group I can display the common names as well. please assist


texts: Finding suitable font size for map scale in ArcGIS Desktop?

texts:
Is there a website that suggests suitable font sizes for each scale level in the layout?

I am currently using ArcGIS 9.3 but pretty soon will upgrade to 10.


texts: Create a script to add delimited text layer in QGIS

texts:
I'm trying to use QgsVectorLayer() function to add delimited text layer in QGIS. So I have a txt file (T_1999_1_2.txt) in path: "D:\PATRICIA\DOCUMENTOS\ESTACOES METEOROLOGICAS\DADOS METEOROLOGICOS\T" and I want to create a shapefile named: T_1999_1_2.shp. The txt file structure is:

  par;num;ano;mes;dia;D;btq
  3;705;1999;1;2;-999,9;S;40,43885833;-8,43994167
  3;619;1999;1;2;25,0;A;41,70972;-8,02699
  3;718;1999;1;2;104,0;A;39,78055278;-8,82096667
  3;766;1999;1;2;28,0;A;38,76620278;-9,12749444
  3;560;1999;1;2;333,0;A;40,71492778;-7,89591667
  3;669;1999;1;2;-999,9;S;39,83950000;-7,47866944
  3;848;1999;1;2;-999,9;S;38,48486667;-7,47291667
  3;555;1999;1;2;8,0;A;41,70655556;-8,80210833

I have tried to run the script but appears the same window error ("Error in file. File cannot be opened or delimiter parameters are not valid") My sript is the following one:

uri = "D:\PATRICIA\DOCUMENTOS\ESTACOES METEOROLOGICAS\DADOS METEOROLOGICOS\T\T_1999_1_2.txt?delimiter=%s&xField=%s&yField=%s" % (";", "x", "y")
vlayer = QgsVectorLayer(uri, "D:\PATRICIA\IG\QGIS\FWI\T\T_1999_1_2.shp", "delimitedtext")

I don't understand, why it doesn't work.


texts: How are section of roads uniquely identified (in both directions) by Google Maps?

texts:
I am wondering whether roads around the world are uniquely identified by a certain value in Google Maps? If so, are they mapped uniquely in terms of their direction? i.e is route A to B is identified differently from route B to A?

In general, I am looking at how the backend database would be structured in terms of uniquely identifying a road, and also if they split it into the two different directions.

Since I don't really know most of the maps, I guess to narrow down the scope any example would work else just take it that I am asking about how Google Maps does it!


texts: Changing direction units to km in Navigator for ArcGIS

texts:
I created a road network using ArcMap 10.6 in Projected Coordinate System (PDO_1993_UTM_Zone_40N - meters). I have one travel mode in this network(car). there is no restriction attribute in my road network and i am using length(meters) as default impedance and i have travel time as well. I have defined the driving directions (km) as follows.

Directions Configured: Yes
General Directions:
  Display Length Units: Kilometers
  Length Attribute: Length
  Time Attribute: TimeInMin
Source Directions:
  Centreline:
    Street Name Fields:
      Primary:
        Name: NAME
        Suffix Type: TYPE

Using Make a Map with your own data tutorial I did the following steps:

I used ArcGIS Pro 2.4 to created a custom basemap
Added the road network layers to this basemap.
Created a new map and added this basemap to it.
Executed the "Create mobile map package" tool to create a package.
Executed "Share Package" to share it on ArcGIS Online.

Now the issue is when i am using it in Navigator for ArcGIS i am getting a route in miles and i am not finding any option to change it to km.

My main concern is when i create and share the mobile map package i want to use km to be default unit for route.


texts: Looping through 16 million records using ArcPy?

texts:
I have a table with 8 columns and ~16.7 million records. I need to run a set of if-else equations on the columns. I've written a script using UpdateCursor module, but after a few million records it runs out of memory. I was wondering if there is a better way to process these 16.7 million records.

import arcpy

arcpy.TableToTable_conversion("combine_2013", "D:/mosaic.gdb", "combo_table")

c_table = "D:/mosaic.gdb/combo_table"

fields = ['dev_agg', 'herb_agg','forest_agg','wat_agg', 'cate_2']

start_time = time.time()
print "Script Started"
with arcpy.da.UpdateCursor(c_table, fields) as cursor:
    for row in cursor:
        # row's 0,1,2,3,4 = dev, herb, forest, water, category
        #classficiation water = 1; herb = 2; dev = 3; forest = 4
        if (row[3] >= 0 and row[3] > row[2]):
            row[4] = 1
        elif (row[2] >= 0 and row[2] > row[3]):
            row[4] = 4
        elif (row[1] > 180):
            row[4] = 2
        elif (row[0] > 1):
            row[4] = 3
        cursor.updateRow(row)
end_time = time.time() - start_time
print "Script Complete - " +  str(end_time) + " seconds"

UPDATE #1

I ran the same script on a computer with 40 gb RAM (the original computer had only 12 gb RAM). It successfully completed after ~16 hours. I feel that 16 hours is too long, but I've never worked with such large dataset so I don't know what to expect. The only new addition to this script is arcpy.env.parallelProcessingFactor = "100%". I'm trying two suggested methods (1) doing 1 million records in batches and (2) using SearchCursor and writing outputs to csv. I will report on the progress shortly.

UPDATE #2

The SearchCursor and CSV update worked brilliantly! I do not have the precise run times, I'll update the post when I'm in office tomorrow but I would say the approximate run time is ~5-6 minutes which is pretty impressive. I was not expecting it. I'm sharing my unpolished code any comments and improvements are welcomed:

import arcpy, csv, time
from arcpy import env

arcpy.env.parallelProcessingFactor = "100%"

arcpy.TableToTable_conversion("D:/mosaic.gdb/combine_2013", "D:/mosaic.gdb", "combo_table")
arcpy.AddField_management("D:/mosaic.gdb/combo_table","category","SHORT")

# Table
c_table = "D:/mosaic.gdb/combo_table"
fields = ['wat_agg', 'dev_agg', 'herb_agg','forest_agg','category', 'OBJECTID']

# CSV
c_csv = open("D:/combine.csv", "w")
c_writer = csv.writer(c_csv, delimiter= ';',lineterminator='\n')
c_writer.writerow (['OID', 'CATEGORY'])
c_reader = csv.reader(c_csv)

start_time = time.time()
with arcpy.da.SearchCursor(c_table, fields) as cursor:
    for row in cursor:
        #skip file headers
        if c_reader.line_num == 1:
            continue
        # row's 0,1,2,3,4,5 = water, dev, herb, forest, category, oid
        #classficiation water = 1; dev = 2; herb = 3; ; forest = 4
        if (row[0] >= 0 and row[0] > row[3]):
            c_writer.writerow([row[5], 1])
        elif (row[1] > 1):
            c_writer.writerow([row[5], 2])
        elif (row[2] > 180):
            c_writer.writerow([row[5], 3])
        elif (row[3] >= 0 and row[3] > row[0]):
            c_writer.writerow([row[5], 4])

c_csv.close()
end_time =  time.time() - start_time
print str(end_time) + " - Seconds"

UPDATE #3
Final update. The total run time for the script is ~ 199.6 seconds / 3.2 minutes.


texts:
I am working on a residential real estate database in MS ACCESS.  I currently have residential property’s tagged by parcel ID’s.  I am trying to run a query to pull in multiple sales data from various tables.  I need to somehow map in street mapping to only pull data for all property’s within a mile radius.  Can anyone explain to me how I can find this mapping?   I can program this in access – I just can’t locate the mapping for this project.
Thanks,
Derrick

PS - location - US, Indiana


texts:
I have downloaded Sentinel 2 and Aster images with the SCP in QGIS. When I try to build a raster or create a bandset using different band combinations, the images turn out black. When I look at the single bands, it shows me the correct image, but the moment I combine the different bandwaves, I get a black image. What am I doing wrong? Can someone please give me advice?
I did the histogram cut stretch function, but this did not work.


texts: header information is missing from image, after writing using opencv

texts:
I am doing complement of an image using open-cv python, while writing the image, open-cv is not writing the header information in the complemented image. my image data are in TIFF format.

i am doing like this...

import cv2
img = cv2.imread("image.tiff", -1)
for i in xrange(0, len(img)):
    for j in xrange(0, len(i)):
      img[i][j] = 65535-img[i][j]
cv2.imwrite("complemented.tif", img)


texts: When computing flow direction, the D8 algorithm is used - are there major differences if others are used?

texts:
I am a student working with Hydrological models.  Within these models we can route water across a DEM (within a delineated watershed or just within an un-delineated DEM) using many different flow routing algorithms: D8, MFD (sometimes FD8), Rho8, DEMON, D-Inf, and so on, with different variations within even the ones I listed possible (original MFD vs. MFD with a power term).

But, when working within a delineated watershed, a flow routing algorithm (primarily the D8) has already been applied to get the flow directions grid - typically the first step in a delineation.

Has there been any research looking at the difference between delineations performed on flow directions grids created using the other algorithms?  I know since these other algorithms allow multiple flow directions, that it might be difficult to visualize this data (unless there is some compounding into one value) - and that that alone might make it impossible.

It doesn't seem like the delineations would be any different - if point y is in the watershed of pour point x, water will flow from point y to pour point x, whichever way it may get there.  The flow may be split into an infinite number of paths (infinite dispersion), but they will still all end up at point x, right?  So even though D8 might not be perfect for predicting actual flow paths, it should produce the same delineation as any other, right?


texts: Create a Non-North/Rotated-Vertically User Coordinate System?

texts:
I need to plot 3d PostGIS data into vertical planes for symbolization and drawing vertical cross-section interpretations for mineral exploration models. I have written a short plpgsql function that takes the PostGIS 3D geometry, moves an xy-coordinate to (0,0,0) then rotates it on the z-axis based upon a given rotation(degrees) and then rotates it vertically along the x-axis.  This new geometry is then filtered/cropped by the new/rotated z-coordinates so just the right plane of data can be viewed.  It works and functions well with QGIS but the downside is that it is not a very update able view for any points/lines/polygons that I draw for any particular cross-section interpretation.

Is it possible to write a proj4/proj6 project that rotates the data on the x-axis?  It would have to be version of the Oblique Mercator to provide the horizontal/z-axis rotation?

I know I am asking for a lot but vertical cross-section interpretation is crucial in the mineral/petroleum industry.

    CREATE OR REPLACE FUNCTION public.xsec_vert_rotate(geometry, double precision, double precision, double precision)
RETURNS geometry AS
$BODY$
    DECLARE
           orig_geom alias for $1;
           x1 alias for $2;
           y1 alias for $3;
           rot_deg alias for $4;
           rot_geom geometry;
    BEGIN
           rot_geom := ST_RotateX(ST_RotateZ(ST_Translate(orig_geom,(0.0::double precision - x1),(0.0::double precision -     y1), >0),radians(rot_deg)),-pi()/2);

           RETURN rot_geom;
    END;
$BODY$
    LANGUAGE plpgsql VOLATILE
    COST 100;
COMMENT ON FUNCTION public.xsec_vert_rotate(geometry, double precision, double precision, double precision) IS 'Cross-Section (Vertical) Rotate: geometry3d, x1/y1 Point Coordinates to move to (0,0,0), Rotation_degress along z-axis to orient data to xy grid system, sets SRID to Zero';
SELECT l.sampleid,
    l.holeid,
    l.sampfrom,
    l.sampto,
    st_z(xsec_vert_rotate(st_transform(st_startpoint(l.bhline_geom), 3158), 375800, 6307000, -65)) AS z_xsec,
    xsec_vert_rotate(st_transform(l.bhline_geom, 3158), 375800, 6307000, 65) AS bhline_geom
   FROM public.lab
  WHERE st_z(xsec_vert_rotate(st_transform(st_startpoint(l.bhline_geom), 3158), 375800, 6307000, -65)) between -50 and 50;


texts: Warning inconsistent extent error for shapefile with .prj?

texts:
I've got a shapefile which already has a projection, and whose GCS is WGS 1984. When I add it to a blank map document I get:

  Warning! Inconsistent extent! error.

The extent of the shapefile is bizarre. How can I clear the .prj file and correct the extent? At this time I don't need actually need it in a specific projection. Here's the extent:

Top: 1647192.701300 dd

Bottom: 388436.828700 dd

Left: -35734.283900 dd

Right: 940661.882900 dd


texts: How do you change units from meters to kilometres?

texts:
Possible Duplicate:
  Using ‘gdalwarp’ in Python script for Grass

I wonder how i could fit my cells to 10 km... Therefore is in meters... and the resolution is different (0.16164). I need this for create a buffer zone. But it is in meters. I will like to learn solving. If I made a project like WS84, it is the same. What should i do? My map is probabilities.

Thanks

g.region -p
projection: 0 (x,y)

zone:       0

north:      70.38082

south:      29.80918

west:       -9.52082

east:       37.35478

nsres:      0.16164

ewres:      0.16164

rows:       251

cols:       290

cells:      72790


texts: Creating new vector layer that fills space between other polygons in QGIS

texts:
I inherited a clipped (500m radius) land cover map containing vector layers of different land cover categories such as buildings, grass, roads, etc from a colleague.

The empty (white) space in the map is actually paved surfaces in reality but is it possible to fill that "empty space" with another vector layer that creates a new landcover polygon class automatically? I'm under the impression that it is possible to create another layer with the symmetrical difference tool but I think I need two input vector files - the one with all the current layers and another layer. Am I supposed to create another vector layer that covers the exact same area (i.e. 500 m radius about a certain point?) and how do I go about doing that?
To summarize, this is what I did to get the difference layer:

Add a vector layer with center-point to existing maps.
Create a buffer layer using this centerpoint. Vector > Geoprocessing tools > Buffer. Use the new vector layer with center-point as input layer.
With the new buffer layer, you can now use the Difference Tool (Vector > Geoprocessing > Difference). Use the buffer layer as your input layer, and the other layer being your main map.
You will have to edit the attributes of the new layer so it corresponds to the same type of attributes in your base map.
Merge the layers using the mmqgis tool (mmqgis > transfer > merge layers), and create a new shapefile with all the layers merged together.


texts: Reproducing Paper Results for Land Cover Mapping of Satellite Imagery?

texts:
I am interested in training a neural network to perform pixel-based classification for land cover mapping of satellite data. In order to learn more about this interesting field I am trying to reproduce the results from this paper:

ASSESSMENT OF CLASSIFICATION ACCURACIES OF SENTINEL-2 AND LANDSAT-8 DATA FOR LAND COVER / USE MAPPING

I currently have 30 satellite images of my area of interest per year for a total of 5 years, and am in the process of creating the training data set.

  Pixel-based classification: Classification is done on a per pixel level, using only the spectral information available for that individual pixel (i.e. values of pixels within the locality are ignored). In this sense each pixel would represent a training example for a classification algorithm, and this training example would be in the form of an n-dimensional vector, where n was the number of spectral bands in the image data. Accordingly the trained classification algorithm would output a class prediction for each individual pixel in an image.

While understanding that the pixel is the fundamental unit of a satellite image and so it is natural to analyze spectral information contained within, I do not understand how this is used in order to process a whole region which may contain 50,000 pixels.

The authors use a 30m data set. This is represents 3 pixels from Sentinel satellite imagery.

My question is: If I train a neural network on a 30m data set to classify between vegetation and buildings, how would it be possible to test the network to classify an image containing, for example 10,000 pixels?

If I train a CNN on images containing only 3 pixels how can I test effectively on larger images such as the one below?

This is an example of a result I would like to replicate, however am unable to understand this concept which is vital for analyzing spectral information:


texts: How to prevent ArcGIS JS map from failing when using a Bootstrap "hidden" class?

texts:
I'm trying to put an ArcGIS Server JS API map within a div with a hidden-sm class in Bootstrap 3, so that the map is not shown on small devices.

The map's parent div is correctly hidden when dragging a large window to become narrower on a desktop browser.

But when re-expanding the window such that the map's parent div displays again, the map is no longer functional.

Note that the div is shown but the map is empty, and the zoom buttons are non-functional.

How can I hide the map on a smaller device, without causing it to fail when the browser is expanded again?


texts: Creating a GCP (*.points) file from a Shapefile (points layer)?

texts:
I have created a points layer which represent Ground Control Points (GCP) for georeferencing scanned topographic maps via 'Add Layer > Add Delimited Text Layer', where I saved coordinates in a .csv format that has the lat-lon of each point, stored in DMS coordinate system (I am using DMS as the scanned topo maps use DMS).

I have managed to load the .csv as a Shapefile. Now, I open up the Georeferencer, load the scanned topo map, and click 'File > Load GCP points', however it requires a Ground Control Point file with the *.points extension, and not a Shapefile.

I would prefer not to manually reassign all the control points. How do I convert a Shapefile to a GCP (.points) file? If that is not possible, how do I create a GCP (.points) file? What file formats can I use to create the required GCP (.points) file?

My goal is to be able to load it into the Georeferencer. I am using QGIS 3.2.1.

When I change the extension of the .shp to .points, it becomes a POINTS file (.points), but QGIS will prompt me the following; "Invalid GCP file. File could not be read."


texts: Converting a single shapefile with multiple features into independent GeoJSON files using ogr2ogr

texts:
Let's say there's a single shapefile that contains multiple features. Most of which are polygons.
How can one convert those features into separate independent GeoJSON files?
Can this be done with ogr2ogr?


texts: Spatial segregation from a raster of racial distribution in the USA

texts:
I am trying to measure the degree of racial segregation at some geographic level (let's say the county level, for simplicity) in the USA. I have a shapefile with these geographic units.

To do so, I've downloaded this map in which there are dots of different colors (4 colors) and each color represents a different race. Each dot is a person. I have georeferenced the map so I have a properly georeferenced raster.

What would be the best way of computing a segregation index at the  level of my geographic units? Given that:

-In each pixel, there can be different colored dots, so that the resulting pixel is a mixture of colors. It is then not trivial for me to infer the % of each race in each pixel/grid
- I am not that familiar with research on spatial segregation, so I could use some help on the measure I have to choose.

I can do it using: QGIS, GDAL, Python or R


texts:
I understand the basic concepts of GPS positioning system and the use of trilateration but I don't understand why it takes so much time to get a fix.

On commercial GPS like TomTom it usually takes about one minute and on my phone it can take more than 5 minutes to be located.

Is it because of the computation time or it needs to collect more data before having a good approximation?

Thank you


texts: Lost relation when importing from OSM to PostgreSQL with ogr2ogr

texts:
I am currently trying to import all administrative boundaries (countries, states, cities, etc) from OSM.

I downloaded a small file from Geofabrik (Corsica), in order to run my tests. I understand that I could run the same operation worldwide using planet-latest.osm.

First, I extracted all boundaries using the following command like

osmosis --read-xml file="corse-latest.osm" --tf accept-ways boundary=administrative --used-node --write-xml out1.osm

Then, I import everything into PostgreSQL using

ogr2ogr -f PostgreSQL PG:"[connection details]" out1.osm -lco COLUMN_TYPES=other_tags=hstore

Which fills up the multipolygon table with the list of cities, and that is exactly what I wanted.

However, if I look for "Ajaccio", there is a lot of useful tags associated with the node 51341115, that do not show up in the multipolygon. Instead, they are associated to a point linked to that node ID.

In the XML file, I have been able to link that point to the relation 73283, inside of which are standing all the ways of the city, and the tags that I actually find in the multipolygon table.

The question is then the following: how can I link the multipolygon to the point in order to extract the tags, using solely the database?

I'm not specifically tied to those tools, I'll use anything that works and that spares me writing code.

Many thanks :)


texts: QGIS Advanced Digitizing on drawing new feature with 60 degree angle

texts:
I am trying to draw a 60 degree line from the long line whenever I draw a long line regardless the angle of the long line.
But right now when I am using advanced digitizing tools and 90 degree means straight. I want the 90 degree to be my long line so that I can easily draw short line at 60 degree.
I am trying to figure up a method to do this few days and come out with nothing.


texts: Load frequency histogram in Print composer

texts:
I'm a QGIS-user and I'm looking for a way to add a frequency histogram in the print composer. The aim of it would be to visualize the distribution of values around the defined class borders of a given layer.

Considering QGIS can load such charts for graduated layers (properties > style > histogram), I thought it would be interesting to import them straight into a project's print composer.

Does anyone have ever found a proper way to have this done ? I've come across similar topics, like this one How to create a histogram legend in QGIS print composer? ... But in my case the main interest would be to have this fixed quickly and automatically, as this specific kind of chart can already be generated by the program.


texts: Finding shortest routes from one point to all others?

texts:
I have lat and long information for 234 locations and I have shapefile of the route network.

I am trying to obtain shortest path from each location to all other locations. For instance for location 1, I would like to find shortest path from location 1 to 2, location 1 to 3, location 1 to 4  ....... location 1 to 234.

I tried to do this in ArcMap with OD cost matrix but the output was a straight line . I would like to obtain routes like I obtain by using "new route analysis".

Is there a way to obtain routes instead of straight line in ArcMap, or by using python?


texts: PgRouting - How to clip links when reaching max costs?

texts:
I have a polyline shapefile representing a road network and a second shapefile containing points. I would like to use PostGIS (presumably PgRouting) to identify sub-networks or service areas radiating from these points.

Essentially, I am hoping to ask the question, "Starting from point X, how far could I walk in any given direction, given a total travel budget of 1 km, following the road network?" The result would be a set of clipped polylines representing the total range of travel possibility, given a 1 km travel budget.

For reference, this GRASS analysis appears to be exactly what I want to do (except I want to do this in PostGIS): http://www.gdf-hannover.de/lit_html/grass60_v1.2_en/node57.html#sec:optalloc

This next example appears to be almost what I want to do, except it seems to answer the question "which nodes could I travel to given a travel budget of X distance?" http://underdark.wordpress.com/2011/02/12/drive-time-isochrones/

The second is not quite the answer I'm looking for, as I want the polylines clipped to my travel distance--I don't care if I make it all the way to a node.


texts: How to catch when layers visibility change in ArcGIS?

texts: I'm trying to regenerate a Legend in a Pagelayout whenever the layers visibility are changed in the TocControl. What is the best way to catch when a layer visibility is changed?


texts: How to add a save and load button in a Python-addin ArcGIS?

texts:
I would like to generate a tool in Python-Addin ArcGIS as the image attached. I also like to generate the save and load button in the tool. My purpose is that when the user chooses the range parameters and then clicks the save button, those parameters can be saved and next time if the user would like to reuse the parameter, they can load it. This is my python script but I do not know which types of python-addin should be put in my script
    import arcpy, arcgisscripting,math
    import pythonaddins
    import  pandas as pd
    import os, glob

    from arcpy import env
    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension("Spatial")
    from arcpy.sa import *

    input=arcpy.GetParameterAsText(0)
    range_1=arcpy.GetParameterAsText(1)
    range_2=arcpy.GetParameterAsText(2)
    output=arcpy.GetParameterAsText(3)

    if (range_2==''):
       outFocalStat = FocalStatistics(input, NbrAnnulus(range_1, range_1, "CELL"),
                           "MINORITY", "NODATA")
    else:
       outFocalStat= FocalStatistics(input, NbrAnnulus(range_2, range_2, "CELL"),
                           "MINORITY", "NODATA")

    outFocalSat.save(output)


texts: Adding time dimension to coverage in geoserver

texts:
I'd like to add a time dimension to my coverage in GeoServer. I've seen quite a few questions on this in this stack exchange but the solutions don't really work for me.

I need to do this via the REST API, without the ImageMosaic plugin and without changing the file names.

I however first tried with the GUI. I am following the instructions in

http://docs.geoserver.org/stable/en/user/data/webadmin/layers.html#data-webadmin-layers-edit-dimensions

I get immediately stuck because it says there is no time information available.

I looked into the other parts of the GUI to see where I could specify the time for my layer. My best bet was the metadata section but it seems this part is irrelevant.

So, my question is, how do I specify the time for my layer (given that I don't want to use image mosaic nor can I change the name of the file).


texts: Open source software for car navigation on iOS platform

texts:
I'm trying to make a car offline navigation app for iOS. Now I analyze OpenStreetMap features what to use. I would like to know if there is open source solution for the following parts of app:

best storing of geodata - do you know Spatialite? Or use standard file system?
rendering OSM (whatever format/own?) - vector rendering, not tile's system
routing OSM, finding shortest/fastest path

I did some experiments but if you have any experiences, i will be glad if you help me with start.


texts: Python: Arcpy: Adding data to a shapefile

texts:
Complete arcpy beginner here...

I have a shapefile with the fields CNTRY_NAME and GDP_HEAD.

I have a list in Python gdppop2002 = [['Country1, 234.1232341], [Country2, 23.2341234],...etc...]

I'm trying to write each float value in the list to the corresponding country row in the shapefile's attribute table. I need to check if each Country in gdppop2002 matches a CNTRY_NAME, and then write the float value in the list to GDP_HEAP. From a couple of hours of googling this appears to require creating an insert cursor, but I'm a little lost. Would anyone be able to point me in the right direction?


texts: ArcGIS find all polygons that intersect with more that one polygon from another layer

texts:
I have two layers with polygons. I want to find those polygons from first layer that intersect with more that 1  polygon from second layer and find those intersected polygons ID numbers: 1 from first layer and all from second layer, that intersect with it. I do not need polygons from first layer that have no intersection or intersect only with 1 polygon from second layer.

Is it possible by using ArcGIS to achieve this? Maybe somehow by selection only needed polygons or by joining all intersected polygon numbers from second layer to the first layer properties.


texts: I have a bit of confusion on the pre-processing of remote sensed images. Do  atmospheric correction methods such as Dark Object Subtraction convert the DN into reflectance values? If not, it is considered an acceptable pre-processing to perform only a Dark Object Subtraction, or before is always required a radiometric calibration?


texts: Troubles translating own OSM server tile to English (ubuntu 18.04)

texts:
I have created a OSM tile server using the https://switch2osm.github.io/ tutorial for ubuntu 18.04 with great success but no matter what I cannot display the map labels in English only, I have tried many solutions from the Open Street Map help forum but none seem to work for me.
I tried adding the line node,way name:en text linear to the default.style in osm2pgsql so the field name:en gets included in the postgres database, but the name:en field never appears on the tables even re-building osm2pgsql, dropping the current table and re loading the map data into db did not help(I am using https://download.geofabrik.de/asia/azerbaijan-latest.osm.pbf from the switch2osm guide).
Maybe I am just being clumpsy and missing something out, if any person with related experience can guide me I would be really grateful.


texts:
Is there an easy way to have the legend linked to the size of the box it is drawn in.

So instead of running down the side of the page, it will be across the bottom in columns?

Does that make sense... I don't want it to run as a line across the bottom either but define my box with lets say 3 items in each column?


texts: How do I convert a QGIS 1.8 plugin to QGIS2.0?

texts:
I have developed a plugin under QGIS1.8
and now I want to make it work on QGIS2.0

I have chage qgisMinimumVersion into 2.0(originally 1.8)

but QGIS2.0 still give me the error say

This plugin is incompatible with this version of QGIS
Plugin designed for QGIS 1.x

how to fix it??


texts:
I've got GDAL installed from the latest 1.11.2 source, but when I try to make by C++ executable I get undefined references to "GDALAllRegister" and basically all the GDAL functions I call.  The same cpp file compiles fine on my Mac, but I can't seem to get the include directories correct.

In my CMakeLists.txt I have:

find_package( GDAL )
include_directories( ${GDAL_INCLUDE_DIRS} )

And I've made sure GDAL_INCLUDE_DIRS points to /usr/local/lib where the gdal libraries are located.


texts: Can I clip a georeferenced tiff to a shapefile?

texts:
If I import a tiff file of a zoning map for a city and then georeference it to a shapefile of the city boundary, is it possible to clip the georeferenced tiff file to the city boundary shapefile?

I have tried to use the extract by mask tool to do this, but even after making sure the tiff and shapefiles were in exactly the same projection, the tool would not work.


texts: R statspat package (or similar pkg) and convert back to spatial objects

texts:
I've been testing the waters and doing some basic overview on the statspat package.  This is a complex package and would likely require some serious effort to jump in and then get familiar with the package.
This package seems to have some great features - including drawing a radial line from a given point at some defined angle. This feature saves the effort from doing the math to accomplish the same task.  And it accomplishes the task with a simple line of code that is intuitive and very readable.
plot(0, xlim=c(0,5), ylim=c(0,5), main=" radial & sp ", xlab="", ylab="", type="n")

    draw.radial.line(0, 2, center=c(3,3), angle=90)         # create spatstat object

# create spatial linestring object and combine with the plot:

    p1 <- c(3,2)
    p2 <- c(3,4)

    mypt <- matrix(data = c(p1, p2), nrow = 2, ncol = 2, byrow = T)
    mypt.sf <- st_linestring(mypt, dim = "XY")
    mypt.sf
    plot(mypt.sf, add = TRUE)

The spatial linesting object and the radial line object both plot nicely together.
But since these two lines are different types of objects, the st_intersects() and other similar functions don't recognize the radial line.  Even though these object plot together, it looks like these objects won't play together well.
The spatstat package has many features that convert spatial objects to spatstat objects, but there doesn't seem to be any similar features for converting spatstat objects back to sp objects.  I looked around and haven't found any way to convert spatstat objects to sp objects.
So this question boils down to asking whether there is a way to convert a spatstat objects to spatial or simple feature object?  Or is there another package that provides some simple functions similar to the spatstat, draw.radial.line() function ? And if so, do these objects convert back and forth and play well with sp or sf objects?
Link to plot link


texts: Using database tables with bigint columns in ArcMap query layer?

texts:
I am very new to GIS and am using the ArcMap 10 program.

Currently, I am trying to set up a query layer from a remote database. I am able to connect to the database, but I am having trouble making the query because the table I want to use has a "bigint" column. From what I have read, ArcMap does not support this data type directly, though some sources suggest that there is a workaround by manually writing a query in the "name" and "Query" fields instead of using the "List of tables" and "columns" selection menus. I've tried this "workaround", and have not succeeded.

Is there any way that I can use my database in a query layer without modifying the database at all? Or can I connect the database to ArcMap in another way?

I am using a PostgreSQL DBMS.


texts: I can´t open a GeoTiff layer in QGis Lyon

texts:
I get a problem.
I have been trying to open a Geotiff layer and Qgis turn off automatically. Then a report appear (bellow).
I hope someone can help me with it.

Dyld Error Message:
Symbol not found: _GTIFAllocDefn
Referenced from: /Library/Frameworks/GDAL.framework/Versions/1.11/GDAL
Expected in: /Library/Frameworks/UnixImageIO.framework/Versions/E/UnixImageIO


texts:
Why can't I add a base map to my data frame in ArcMap?

I am using ArcGIS 10.2 for Desktop and the options under Add Data for Add Basemap and Add Data From ArcGIS Online are both disabled (grayed out).


texts: ST_LocateAlong interpolation: difference between PostGIS and Spatialite

texts:
In a previous question, I discovered that ST_LocateAlong along a LineStringM is not implemented the same way in Spatialite than in PostGIS.

While PostGIS interpolates between measures, this is not the case in Spatialite. In Spatialite, only vertices which exactly match the m-values are returned.

Example:

SELECT
    ST_AsText(ST_LocateAlong(ST_GeomFromText('LINESTRINGM(0 0 1, 0 1 2, 1 1.5 3, 1 2 3.5)'), 2));

PostGIS returns MULTIPOINT M (0 1 2)
Spatialite returns POINT M(0 1 2)

Both implementations return (almost) the same result, except that PostGIS is MULTIPOINT M instead of POINT M.

When the measure is not on a vertex, PostGIS interpolates it while Spatialite returns NULL:

SELECT
    ST_AsText(ST_LocateAlong(ST_GeomFromText('LINESTRINGM(0 0 1, 0 1 2, 1 1 3, 1 2 3.5)'), 2.5));

PostGIS returns MULTIPOINT M (0.5 1 2.5)
Spatialite returns NULL

Is there an (easy) way to emulate PostGIS's behaviour in Spatialite?


texts: Clipping Rasters with a mask layer error

texts:
I've got a NDVI layer that I'm trying to clip with a vector file but when I try and use GDAL Tools clipper I get an error.  That shows up like this:

What am I doing wrong?  This is before I even pick the mask layer!

Traceback (most recent call last):
  File "C:/PROGRA~1/QUANTU~1/apps/qgis/./python/plugins\GdalTools\tools\widgetPluginBase.py", line 28, in someValueChanged
    self.emit(SIGNAL("valuesChanged(const QStringList &)"), self.getArguments())
  File "C:/PROGRA~1/QUANTU~1/apps/qgis/./python/plugins\GdalTools\tools\doClipper.py", line 105, in getArguments
    return self.getArgsModeMask()
  File "C:/PROGRA~1/QUANTU~1/apps/qgis/./python/plugins\GdalTools\tools\doClipper.py", line 143, in getArgsModeMask
    if self.alphaBandCheck.isChecked():
TypeError: 'sip.methoddescriptor' object is not callable

Python version:
2.5.2 (r252:60911, Feb 21 2008, 13:11:45) [MSC v.1310 32 bit (Intel)]

QGIS version:
1.7.3-Wroclaw Wroclaw, 00624b3

Python path: ['C:\Users\dean.walker\.qgis\python\plugins\GeoCoding', 'C:\Users\dean.walker\.qgis\python\plugins\GeoCoding', 'C:/PROGRA~1/QUANTU~1/apps/qgis/./python', 'C:/Users/dean.walker/.qgis/python', 'C:/Users/dean.walker/.qgis/python/plugins', 'C:/PROGRA~1/QUANTU~1/apps/qgis/./python/plugins', 'C:\PROGRA~1\QUANTU~1\bin\python25.zip', 'C:\PROGRA~1\QUANTU~1\apps\Python25\DLLs', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\plat-win', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\lib-tk', 'C:\PROGRA~1\QUANTU~1\apps\qgis\bin', 'C:\PROGRA~1\QUANTU~1\apps\Python25', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\site-packages', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\site-packages\win32', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\site-packages\win32\lib', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QUANTU~1\apps\Python25\lib\site-packages\wx-2.8-msw-unicode', '/usr/share/qgis/python', 'C:\PROGRA~1\QUANTU~1\apps\qgis\python\plugins\fTools\tools']


texts: Reading shapefiles using for loop: Error in readOGR() : length(dsn) == 1L is not TRUE

texts:
I'm using a for loop to read in multiple shapefiles within a directory and manipulate them individually. Presently I cannot do this because the for loop results in the error: "Error in readOGR(shp_paste) : length(dsn) == 1L is not TRUE". Not a single shapefile is read in; this error occurs immediately. I do not know what this means, and can only find a similar post with this error for writeOGR (link here). When I try to load the shapefiles individually into R Studio, I'm successful in doing so (though I haven't tried individually for all of them, there are too many). What could be causing this error? It's difficult to go through the files in the directory individually to identify potential problems because there are over 10,000.
Below is my code:

    library(rgdal)
    input_path<- "/Volumes/Directory/data_folder/"
    files<- list.files(input_path, pattern = "[.]shp$")

    for(f in files){
    ifile<- list.files(input_path, f)
    shp_paste<- paste(input_path, ifile, sep = "")
    shp1<- readOGR(shape_paste)
    print("success!")
    }

Running this code produces the error above.


texts: ArcGIS 10: Adjust custom grid intervals based on map scale?

texts: I am preparing a pmf which includes a custom grid overlay. To clarify: a custom grid overlay is marginalia, not a feature class with geometry (so a scale dependecy inside a layer won't work).  The final user, in ArcReader, will interact with this map at any scale, zooming in and out as they please.  So far, the grid intervals stay the same, regardless of scale, which results in too many grid values at small scales and too few grid values at large scales.  Any suggestions?  Thank you!


texts:
I'm trying to create an Add-in in ArcGIS and used Python Add-in Wizard:
http://www.arcgis.com/home/item.html?id=5f3aefe77f6b4f61ad3e4c62f30bff3b

Everything went fine until i came to this step:

when i click on 'Install Add-In' I get the following Error:

After some digging I came to the conclusion that the wizard creates Addins for 10.1 or above and I am using 10.0. Is there a way to copy-paste the .py file and make it work?

I tried to create a folder at:
C:\Users\my_user\OneDrive\Documentos\ArcGIS\AddIns\Desktop10.0\my_plugin
with the .py file but it doesnt apear in the 'Add-In Manager'.

I made it work following this:
ftp://ftp.dnr.state.mn.us/pub/gisftp/2014_GISLIS/Intro_To_Addins/Pre-Class%20Software%20Setup%20Instructions_AddIns2014.pdf

where instead of VS 2010 Express i used Visual Basic 2008 SP1


texts: Easier way to get Google Maps Directions data into smartphone?

texts:
I often use Google Maps Directions to draw routes before taking my bike, but it's a bit of a hassle to turn its data into a KML file and copy it into my smartphone:

Go to Google Maps Directions https://www.google.fr/maps/dir/
Draw route
Get and copy URL through the Menu button
Go to GPSVisualizer, paste URL, and get KML file
Copy KML file into smartphone (through either Windows Explorer or iTunes).

Is there an easier way, especially for non-geeks among my friends?


texts:
I want to draw a rectangle. The coordinates are somewhere in northern Germany. But the rectangle shows up in the Near East? What am I doing wrong?

       var epsg900913 = new OpenLayers.Projection('EPSG:900913');
       var epsg3035   = new OpenLayers.Projection('EPSG:3035');

       var p11 = new OpenLayers.Geometry.Point(4281000, 3502000).transform(epsg3035, epsg900913);
       var p12 = new OpenLayers.Geometry.Point(4595000, 3502000).transform(epsg3035, epsg900913);
       var p13 = new OpenLayers.Geometry.Point(4595000, 4002000).transform(epsg3035, epsg900913);
       var p14 = new OpenLayers.Geometry.Point(4281000, 4002000).transform(epsg3035, epsg900913);
       var p15 = new OpenLayers.Geometry.Point(4281000, 3502000).transform(epsg3035, epsg900913);


texts: How to get work time-series with AM and PM data?

texts:
I read tutorial and get all work with data like 20170526, but how to get it's work if I have AM and PM data for specifieded date?

As I understand from docs date placing and shp to date field, and then specifying in URL in &time=2017-05-26 format. But as far as I know date do not support AM and PM dates.

Is there any solutions?


texts: Poor resolution / blurry geotiff raster layer on Geoserver with openlayers gui using wms

texts:
I am having a problem with blurry geotiff raster layers served by geoserver and viewed in openlayers browser gui using wms.

I am hoping someone can point me in the direction of a fix.

Openlayers is configured to use the same SRS as the geotiff for raster, ESPG: 4326
SRS for layer confirmed to be projecting at 4326 in geoserver
I have tried the following WMS settings:
disabling jpg & png compression
applying superoverlays
with and without antialiasing
all interpolation options: nearest neighbor, bicubic, bilinear.

The result is the same. blurry / low resolution image is present in openlayers web gui as well as geoserver openlayers preview.

one odd thing, the layer's resolution is correct on the last two step/zoom levels. IT looks the same in native preview and openlayers, But at all other zoom levels there is blur.

Please see screenshots for reference. The image on the left is gotiff in native image viewer. openlayers gui is on the right.


texts: Dissolve a polyline feature class so that touching features dissolve into a single feature

texts:
I'm trying to dissolve a polyline feature class so that spatially intersecting features become a single feature. When I run the Dissolve tool, setting Create multipart features and Unsplit lines both to false, I get unexpected results. In the example below, two line features I'd hoped to dissolve into one became three line features:

Is there a way to dissolve a line feature class so that any touching features dissolve into one?


texts: Use of consistent symbology across multiple stretched rasters in ArcGIS Pro

texts:
I have multiple rasters that represent the same thing (NDVI) at different times. Each of the rasters has a different range of values. I want to symbolise them using the same range of values (e.g. from 0 to 1) so that they can be compared. What is the best way of doing this?

This question has been answered for ArcGIS. However, the procedure is different for ArcGIS Pro as the options have been removed / changed.

The image below shows three of the layers with their differing value ranges:


texts: Project tiles from custom imagery in transverse mercator with OpenLayers

texts:
I have the following parameters of a custom imagery and I need to set-up a layer in OpenLayers with which I would project the tiles properly from this imagery.

TransverseMercator
ellipsoidName=WGS84
equatorialRadius=6378137.0
inverseFlattening=298.257223563
falseEasting=500000.0
falseNorthing=0.0
origin=[latitude = 0.0 longitude = -117.0]
scale=0.9996

I apologize that I did not provide you with enough information. I do not understand how projections work and will appreciate any information that could lead me to a solution.

Do not hesitate to tell me what other information would be useful.


texts: Using ArcPy to verify paths of Hyperlinks?

texts:
I've recently come across an issue where some of our hyperlinks aren't working.
I'm familiar with python and have made similar scripts before, but I would first like to ask for your assistance to maybe something already pre-existing or another tool to use.

I would like to cycle through all the hyperlink paths for a certain shapefile and search the C:\Scans\ for the path.
If the path doesn't exist than I want the path to be written to an excel file.

Can anyone help me with this?

Final code so far

import arcpy
import os

fc = r"G:\GIS\Lauren\Shapes\watermains.shp"
fields = ["Truck_Path"]

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


texts: I have some MODIS data processed with MRT tools. They are now in tiff format. Now, I have to complete Density slicing to identify water pixels from the image. The DN higher than 1985 is usually non-innundated area where lower than that value is the inundated area. I completed the task with ILWIS but want to do that with QGIS. Is it possible to do that? Water pixel will be red color and non-inundated area will be green color.


texts: I'm curious what solutions there may be for generating Mapbox vector tiles from SDO_Geometry data? There doesn't appear to be native support as with PostGIS.


texts: Using python to convert string values with Calculate Field in ModelBuilder?

texts:
In ArcGIS 10.5, I am building a model that should convert a field of coded values to another code. This seems like it should be a simple if/then statement but I am not getting it. To simplify the process, I populated the target field with the source field's values and then intend on using python to transpose the values. It should be considered a string, but the values are numerical. The field I am trying to edit is "TP_Condi"

In the "Calculate Field" code block:

def reclass:
if TP_Condi = 3:
     return 1
elif TP_Condi = 2:
     return 8
elif TP_Condi = 1:
     return 3
else:
     return "error"

This is the code I was seeking:


texts: Filling holes in polygons automatically using QGIS

texts:
Is there an option in QGIS to fill holes in polygons automatically by setting parameters e.g. maximum area, distance between borderlines?
I want to delete/fill holes that a very "thin"/small.


texts: Is it possible to create a round map frame in QGIS Map Composer to create circumpolar maps?

texts: Does anybody know if it's possible to create a round map frame in QGIS Map Composer to create circumpolar maps?


texts:
Swath of Sentinel image is 290km. But the images which we download comes only at 190km.

Where can be 290km images be downloaded?


texts: Understanding Inverse Distance Weighting Bulls Eye Effect?

texts:
Is the Bulls Eye Effect generally caused by a sample point with a high value being too far from other sample points?

Or rather a sample point with a higher value than others close to it?

What exactly causes this?


texts: Domain cannot display existing values when editing through the domain

texts:
I have set a series of domains for my Feature Class in ArcMap. The domains work well but it's just when the FC is under edit, it does not display the existing value sitting in the database.

Instead it displays a blank or something similar to "keep existing value".

Is there a way for the domain to display what the existing value is, as well as still have the domain function?

Do I need to first allocate to domain values first in order for the values to be read whilst editing?


texts: Using ST_LINE_INTERPOLATE_POINT in OGR and Python

texts:
I'm interested in finding the midpoint on a linestring using OGR and Python. This can be a line with more segments. Could I use the ST_LINE_INTERPOLATE_POINT in Python using the ExecuteSQL() method?

If so, could someone perhaps help me formulate the correct SQL query?
Or is there an alternative option? I considered using Centroid but the centroid can also lie outside of the geometry.

---- UPDATE ----
I have QGIS installed on Windows and use the OSGeo4W shell. I made the following Python code but unfortunately it doesn't work, geom is seen as a NoneType.

from osgeo import ogr

outdriver=ogr.GetDriverByName('SQLITE')
source=outdriver.CreateDataSource('temp')
tmp=outdriver.Open('temp')
layer = tmp.ExecuteSQL("""SELECT ST_AsText(ST_Line_Interpolate_Point(ST_GeomFromText('LINESTRING ( 90 379, 99 278, 189 224, 297 193, 403 189, 459 202 )'),0.5))""")
geom = layer.GetNextFeature().GetGeometryRef()
print geom.ExportToWkt()


texts: Connecting Python Script tool as Input in ModelBuilder?

texts:
I'm trying to connect my Python Script (it's already in my toolbox and working) with the rest of my Model. So the Script should be my input, but it is not working. I don't know why.

This is how the script looks like:

import arcpy, os
from arcpy import env
env.workspace = arcpy.GetParameter(0)
Topic = arcpy.GetParameterAsText(1)
Alternative = arcpy.GetParameterAsText(2)
featureclasses = arcpy.ListFeatureClasses(feature_dataset='A_'+str(Topic))
Output = arcpy.Union_analysis(featureclasses,r"C:\01_Data\15_Test.gdb\V_"+str(Alternative)+"_"+str(Topic),"NO_FID")
arcpy.SetParameterAsText(3, Output)

The output is set on Type = Derived and Direction = Output.

I tried also this:

import arcpy, os
from arcpy import env
env.workspace = arcpy.GetParameter(0)
Topic = arcpy.GetParameterAsText(1)
Alternative = arcpy.GetParameterAsText(2)
featureclasses = arcpy.ListFeatureClasses(feature_dataset='A_'+str(Topic))
result = arcpy.CreateFeatureclass_management(r"C:\01_Data\15_Test.gdb", ""V_"+str(Alternative)+"_"+str(Topic)")
Output = result.getOutput(0)
arcpy.SetParameterAsText(3,Output)
arcpy.Union_analysis(featureclasses,r"C:\01_Data\15_Test.gdb\V_"+str(Alternative)+"_"+str(Topic),"NO_FID")


texts: Web mapping service with GEBCO bathymetry data

texts:
Is there any open web mapping service provide map tiles with GEBCO bathymetry data. As far I know Only Google uses GEBCO bathymetrydata for ocean floor model. Others uses ocean as plain blue. I need a service provider of web mapping tiles which provide sea floor with GEBCO data.

Thanks


texts: How to remove white parts in ArcMap 10.1?

texts:
This is the part of my map.

On this picture you can see white strips. Some of this white parts are appeared on the intersection of 2 Mosaic Dataset (for eaxample, there is the real map in 1 Mosaic dataset, there is  the empty scan in 2 Mosaic dataset). How to set option "no-data" of Mosaic Layer? I want to display map without white parts.

UPDATE

On this picture you can see different MosaicDatasets. On the other forum people tell me that it is impossible to delete white parts in this case. (because it is not a single mosaic dataset)

Is it true?


texts: Invalid MXD Filename error in one script, but not the other

texts:
I have created 2 geometric trace models for a project, one to pick a spot on the network and run an isolation, and the other to choose from a list of buildings, then then run a trace to find features to isolate that building. Both models use a python script at the end to combine all outputs from multiple group layers into one out feature class for each participating feature class in the network. Also, both models reference the same MXD file in order to used some arcpy.mapping functions.

I published the network iso model on the server and it runs just fine with no errors both on ArcMap and the WebApp I made on ArcGIS online. When i publish the building iso model, it fails with a "invalid MXD Filename" error. The whole output can be seen below.

Both models uses the exact same script, yet it only fails on the building iso model and not the network iso model. I even tried saving a copy of the MXD as a different name, but still received the same error. the part of the script that is failing can be seen below.

The filepath does not have any special characters in it, and again, works just fine in the other model. Has anyone ran into this problem before? I have been looking around blogs and ESRI resource sites but haven't found anything to help. One odd thing i did notice is that the line the log reports the error occurring on is not the same line number as it is in the original script, but i'm not sure if that means anything.

@KHibma provided help. I found out the script is not copying the MXD file to the server, even though it say it is when publishing. Also, it does not seem  to be setting the variables correctly on the server. in the picture below, the script on the left is the network iso, and the right is the building iso.

Why isn't the publishing tools correctly formatting the script from the building iso model?

All the parameters are set the same as in the network iso model.