texts: PyQgis display vertex marker when in editing mode

texts:
Documentation says when issue the following command on a vector layer (polygon), the vertex marker will be displayed in editing mode:

layer.startEditing()

But it is not displaying the vertex markers.

I found that there is a static drawVertexMarker method on QgsVectorLayer class, but do I have to manually call it to draw a marker on each vertex? Also, if I am moving vertex in editing mode, do I have to constantly re-draw all vertices?

Additionally, there is the editBuffer() property that returns the QgsVectorLayerEditBuffer reference, is it possible to display that buffer layer in editing mode for visual aid?


texts: How to count point attributes that fall within a polygon in postgis?

texts:
I have a point shapefile called pointA that has a field called fieldB, and a polygon file called polygonC, and I want to calculate the sum of fieldB for all the points that are within polygonC

How to do this?


texts: sql server geometry find points intersection buffer

texts: How do I find the point of intersection of buffer?


texts: Database Connection error: Invalid database name

texts:
I just created a new database in my Postres called GlenMary and I can't get ArcMap to connect through my already established database connection in ArcMap. I am trying to add some new data to the database but I keep getting an error when I try to connect.

Even though I can connect to the database below it ("beauregard"-highlighted), and the database name validates and populates in the database dropdown, when I click "OK" I get the following error. Any idea how to fix this?


texts: Different results given by PostGIS and MySQL?

texts:
I'm developing a spatial database of administrative areas. I started out with PostGIS, but as MySQL is being used in our other projects, I thought I would at least try to create a spatial-enabled MySQL database.

However, spatial queries give me different results than PostGIS does. Let's say I have a query like this

SELECT nazev FROM orp as o, kraje as k WHERE
  Intersects(o.geom,k.geom) AND Within(Centroid(o.geom),k.geom) AND
  k.nazev = 'Nazev'.

It should return a list of administrative areas within a county with defined name (nazev). Instead, it gives me expected areas + areas that touch boundary of a county, but are not contained by this county. When I run the same query with PostGIS, I get the correct results.

Have anyone ever faced this issue? I am using a national CRS (EPSG:102067) and have tried also with EPSG:4326 with same result. Isn't the implementation of the OGC spec the same as in PostGIS? Thanks for any hint.


texts: Leaflet + GeoServer WFST - Inserting Polygons to PostGIS database

texts:
I have a GeoServer installation on top of a PostGIS database. I am trying to set up a leaflet webmap (using Leaflet WFS-T plugin) so that I can insert, update and delete polygon features. I have managed to get this to work fine for Point features but something seems to be going wrong when inserting polygons.

Below is the request payload to GeoServer:

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0">
<wfs:Insert><schememapper:geoserver-polygons xmlns:schememapper="schememapper">
<schememapper:geom><gml:Polygon xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326" srsDimension="2">
<gml:exterior><gml:LinearRing srsDimension="2"><gml:posList>-1.3735914230346682 52.03958374014823 -1.9696426391601565 50.929466759820585 -3.7992525100708012 51.977294016543425 -1.3735914230346682 52.03958374014823</gml:posList>
</gml:LinearRing></gml:exterior></gml:Polygon></schememapper:geom></schememapper:geoserver-polygons></wfs:Insert></wfs:Transaction>

And this is the response:

<?xml version="1.0" encoding="UTF-8"?><ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://35.242.177.142:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="InvalidParameterValue">
<ows:ExceptionText>Error performing insert: java.lang.String cannot be cast to org.locationtech.jts.geom.Geometry</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

When I look at the geometry type in PostGIS it says it is MultiPolygon so I am conscious that the WFST request is only specifying 'polygon'.

Below is my Leaflet code (I am using Leaflet draw to create the polygons):

var wfstPoly = new L.WFST({
    url: 'http:localhost:8080/geoserver/ows',
//typeNS is the workspace in geoserver
    typeNS: 'schememapper',
    typeName: 'geoserver-polygons',
    crs: L.CRS.EPSG4326,
    geometryField: 'geom',
    style: {
      color: 'blue',
      weight: 2
    }
  }).addTo(map)
    .once('load', function () {
      // map.fitBounds(wfstPoly);
    });

// Drawing
var drawControl = new L.Control.Draw({
    draw:{circle:false, circlemarker:false, rectangle:false,
          },
    edit:{featureGroup: wfstPoly } });
map.addControl(drawControl);

map.on('draw:created', function (e) {
    var layer = e.layer;
    wfstPoly.addLayer(layer)});

map.on('draw:edited', function (e) {
    var layers = e.layers;
    layers.eachLayer( function (layer) {
        wfstPoly.editLayer(layer);
        });
});

// Save button
L.easyButton('fa-save', function () {
         wfstPoly.save();
     }, 'Save changes').addTo(map);

Furthermore, when I use the console to look at the original features pulled from GeoServer they are GeoJSON but when I look at the newly created features they don't seem to be. They just have an array of coordinates.

EDIT: This may have something to do with GML versions. I.E. Geoserver is expecting the XML of the request to be structured differently.


texts:
Is it possible to have ArcGIS 10.0 dekstop log the SQL it sends to the database (direct connect with Oracle 11gr2)?

I am trying to diagnose some performance issues and would like create some scripts to duplicate what this application is sending to the database. However, you cannot obtain those with a SQL trace (complex datatype).

If not - I will capture the network traffic instead, although that less desirable as it is a rather manual process.


texts: Using Google Map API for Business with Openlayers

texts:
Has anyone used the paid version of Google Map API (for business) with OpenLayers? I am currently using the free version and wonder what the difference is, from both the end user and developer perspectives?

I've read the terms on Google site but couldn't find anything specific enough for working with OpenLayers. For example, does it support higher zoom level with clearer satellite imagery?

Thanks!


texts:
I keep getting an insert error from geoserver when trying to save a point, poly etc etc to my postgis table.  any ideas, my code is below

    var Projection = new OpenLayers.Projection("EPSG:4326");

    var map = new OpenLayers.Map("Map", {
        //projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection:Projection,
        controls:[]
       }
    );

    //add google map
    var GoogleMap = new OpenLayers.Layer.Google(
        "GoogleSatellite",{
            type:G_SATELLITE_MAP,
            isBaseLayer:true,
            sphericalMercator:true,
            maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34),
            numZoomLevels: 21
        }
    );

    var NYC = new OpenLayers.Layer.WMS(
      "nyc_buildings - Tiled", "http://soiisdev:8080/geoserver/wms",
        {
            layers: 'cite:nyc_buildings',
            styles: '',
            format: "image/png",
            tiled: 'true',
            transparent: true,
            Projection:Projection
        },
        {
            buffer: 0,
            displayOutsideMaxExtent: true
        }
   );

   var saveStrategy = new OpenLayers.Strategy.Save({auto:true})

    var vectors = new OpenLayers.Layer.Vector("Vectors", {
        strategies: [
            new OpenLayers.Strategy.BBOX(),
            saveStrategy
        ],
        Projection:Projection,
        protocol: new OpenLayers.Protocol.WFS({
            version: "1.1.0",
            url:  "http://soiisdev:8080/geoserver/wfs",
            featureNS : "http://www.opengeospatial.net/cite",
            featureType: "nyc_vectors",
            srsName: "EPSG:4326"
        })
    });

    map.addLayers([GoogleMap, NYC, vectors]);
    //map.zoomToMaxExtent();

    var LonLat = new OpenLayers.LonLat(-73.988113, 40.750898);
    LonLat.transform(Projection, map.getProjectionObject());
    map.setCenter(LonLat, 15);

ive checked the header in fiddler and below is the xml sent to geoserver

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0"  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance xsi:schemaLocation="http://www.opengis.net/wfs  http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
<wfs:Insert>
 <feature:nyc_vectors xmlns:feature="http://www.opengeospatial.net/cite">
   <feature:the_geom>
     <gml:Point xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
       <gml:pos>-73.986374928559 40.742379645006</gml:pos>
     </gml:Point>
   </feature:the_geom>
 </feature:nyc_vectors>


texts: Can you have layer control in Odyssey.js?

texts:
I've been playing around with odyssey.js and made a few map stories with it:

http://www.burdgis.com/blog/archive/map-stories/

What I can't work out is how to toggle layers. I imagine you would need more than one viz.json to do this but only the first one in my code seems to be read and if I add another later on in the article it's ignored.


texts: Summarize the total area of land use polygons within map sheets - how?

texts:
I would imagine that at least some of you have encountered similar problems to mine and I´d like to hear how you solved them.

I have a map sheet index (a grid as a shape file) covering a large area. I also have land use polygons covering roughly the same area. Now I would like to calculate how many hectares of land use polygons fall inside each map sheet (i.e. only the portion of the land use polygons that fall within the map sheet polygon should be taken into account when summarizing "land use area"). The picture below might help to see what I mean, map sheets in blue lines, land use parcels in light green.

Earlier the easiest way to do this for me was to use Hawth´s Tools polygon-in-polygon analysis (http://www.spatialecology.com/htools/polypolyanalysis.php), it did exactly what I wanted. Hawth´s Tools is now replaced by GME and the polygon-in-polygon analysis tool in GME no longer functions the same way as the old one. The new tool only sums up the total area of all land use parcels that intersect each map sheet and writes that as the result. No good for me.

I could use the map sheet lines to cut the land parcel dataset into small bits and then re-calculate the area for all polygons, but that´s not that tempting: I have 26 000 map sheets and 1,2 million land use parcels.

I´ve tried this in ArcGis Basic but no luck there. OpenJump, Qgis or gvSIG do not seem to have a solution either. Do anyone of you?


texts: Matching one layer object color with another layer object color in QGIS

texts:
I'm kind of very beginner in QGIS. I have made one layer of "cells" and visualised them on map with different colors representing cell numbers like this:

The idea is that this layer is like database of the network which is always loaded in my QGIS project and sometimes I need to add some data from these cells on the map, which is loaded in another txt file(s).
I would like the points of this data to represent the same color which is stored in another layer(network database layer). Is there any way how to match these colors?
For example here below you can see the image of the 2 layers loaded together and I added labels to the points. For example how can make points with attribute number 280 to be orange or number 281 to be green?


texts:
I am trying to create a spatial reference using ImportFromESPG and this is what is returned:

AttributeError: 'SpatialReference' object has not attribute 'ImportFromESPG'

Code and error are copied below. I copied this code exactly from the cookbook: https://pcjericks.github.io/py-gdalogr-cookbook/projection.html
so it seems that something must be wrong with my installation?

Windows 10, Python 2.7.8 with IDLE. I get a very similar error when I use python in the command window where the canopy install of python is my default (python version 2.7.11).

Python 2.7.8 comes from the ArcGIS install (10.3.1 with 64bit geoprocessing update).  I installed GDAL from the gisinternals website. Both versions of python on my system are MSC v.1500 64 bit. For Canopy, gdal is installed internally as a package.

Code and error:

from osgeo import osr
sr=osr.SpatialReference()
sr.ImportFromESPG(29613)

Traceback (most recent call last):
  File "", line 1, in
    sr.ImportFromESPG(29613)
AttributeError: 'SpatialReference' object has no attribute 'ImportFromESPG'


texts:
What is a method for importing/exporting DXF files with QGIS?
I work in an engineering field and we use CAD and specialist design programs, and often go through several iterations during a project.
I am having trouble interchanging data between QGIS and DXF formats.  I export several layers from QGIS in DXF format and load them into cad.  Essentially I xref the exported DXF files into CAD, and then draw over the top of them in CAD, as I find it quicker and more accurate.  The files line up OK, but need to be scaled and transformed to match the coordinates in QGIS.  I then detach the xref's, save the new DXF with line work on a single layer.  The DXF file coordinates match the QGIS coordinates.
I can't seem to get the new DXF to import back into QGIS.  I would like to be able to export and import a particular data table several times as the project develops.
Screenshot 1 - I exported the pink polygons and edited them in CAD.  I then imported the DXF file back in but nothing came in - see "entities".

Screenshot 2 - Attribute table of imported layer.


texts: How to delete isolated pixels in ArcGIS 10.1?

texts:
I want to know that it is possible or not to edit the raster pixel in ArcGIS.
I want to delete isolated pixel and I want to delete some pixel by manually.
please tell me how can i do it manually or application of tools.


texts: Pixel Location of Sentinel 1 & Sentinel 2 and SRTM

texts: I am processing integration S1, S2 and SRTM, since their pixel grid does not align with each other I am figuring out how to snap with them all into one grid. Before coming with specific reprojection or resampling techniques I want to know which corner of the pixel is representing ground coordinates for each satellite image (S1,S2, SRTM). I have checked the user guide and technical guide but didn't get the idea yet at the pixel level.


texts: In ArcGIS 10.1, for what tools does Environments actually work?

texts:
As a new user, often I have used "active" (that is, not grayed out and non-functional) Environments when using ArcToolbox Tools and the operation is non-responsive. I wish I documented for my own record what I have seen along the way, but have not, still being somewhat confused as to whether it is because I am doing something wrong. Following are two recent examples.

1) Using a raster, I have seen where

Spatial Analyst > Extract by Polygon > Environments > Raster Analysis
and changing the cell size will not change the cell size;

and

2) using a raster

Data Management > Resample > Environments > Raster Analysis >
Mask will change the cell size but not mask.

Since the combination of operations that are not allowed by ArcGIS do not result in grayed out inactive buttons or a warning message, is there some list that tells what Environment operations can and can't be done with the various ArcToolbox Tools?


texts:
In QGIS 2 the Option for bold Layernames was here:

Go Settings->Options->Canvas & Legend and set 'Bold group names' to true and 'Bold layer names' to false.

This link helps provide context. QGIS 2 bold Layer names?

How is that working for QGIS 3.4 ?
Couldn't find that option somehow.


texts: Combining polylines which share endpoints using ArcGIS Desktop?

texts:
I am working with a shapefile of streets which include links between two intersections represented as two polylines (depending on presence of a median, for instance). I would like to merge such parallel lines with overlapping end points into a single polyline in order to develop a consistent segment-level unit of observation.

I have tried using different tools in ArcGIS Desktop ("dissolve", "unsplit", "merge divided roads", "collapse dual lines to centerlines") and R (gLineMerge) based on other Q&As with similar sounding issues, but not of these tools/packages resolves the issue across the network.

Let's say the input network looks like the image below. The lines that are highlighted are the lines which are divided, but share common begin and end points.

What I would like is for to simplify the network so as to have these divided roads replaced by a single line, as illustrated by the handdrawn coloured lines below.

Any thoughts on what steps could I do in ArcGIS Desktop to automate such process?

The attributes of the lines include "from node" and "to node" columns, but I would like for the geometry to be averaged out ideally (as opposed to removing duplicates).


texts:
I have 2006 land cover GeoTiff from the USGS site (for a county in VA, USA).  The projection is Albers Conical Equal Area Units. The download came with an HTML metadata file that tells me the geo extent in terms of E/W and N/S bounds.  If it matters, I created the dataset using the USGS download tool by selecting a particular county.

I am using GDAL to access the data.

If I take the raster sizes and convert the four corners to WGS84 geo coordinates, I get different geo values for the same pixel value.  In other words, if the x pixel coordinate is 0, then I get a different latitude depending upon whether y is 0 or 4194 (the raster height).  Is this because the rectangular raster image is mapped into a trapezoid on the sphere?

If I take the geo extents from the metadata file and iterate over them, stepping by 1% of the range, converting each geo point to a pixel and line index, I get lots of invalid (negative or too large) index values.  Why would this be?


texts: Intersection and area calculations in geographic vs. projected coordinates

texts:
I think one should do intersection calculations when layers are in the same coordinate system, and one should do area calculations in the appropriate projected coordinate system.

The order matters, but is one preferred? In my system, I get slightly different areas if I

1) Perform intersection of two layers that are in the same geographic coordinate systems, project the resulting intersection to a projected coordinate system, and then calculate areas, or

2) Project two layers to the same projected coordinate system, intersect them, and then calculate areas of the resulting intersection.

Without getting into details of my specific situation, I wonder if one of these approaches is theoretically preferable?  I assume the 2nd would be, but I can't articulate why that would be.

(I assume terms are standard, but this is the intersection I'm looking for: http://wiki.gis.com/wiki/index.php/Intersect)


texts: Changing order of parts in multipart geometry using ArcGIS for Desktop (or ArcPy)?

texts:
Is it possible in ArcGIS to change order of the parts in multipart geometry?

In the context menu i founded "insert before" , "insert after" but it never comes active ... or is it possible via python?


texts: Problem with reprojecting raster from MGI 6 to WGS

texts: When I reproject a raster from Mgi 6 to WGS 84 using Warp in Quantum gis it makes mistake ( about 380m). Is there any way to avoid this?


texts: Sourcing free live traffic layers to use in MXD of ArcGIS Desktop?

texts:
I am looking for free live traffic layers for an ArcGIS Desktop 9.3 MXD.

My searches haven't been that successful.


texts: How to allow the user to place a point as input for a geoprocessing model?

texts: Is there an analogue to Esri's "Feature Set" in QGIS modeller.  I am interested in creating a model that allows the user to place a point feature on a map as an input variable for a geoprocessing model.  I've done this in ModelBuilder, but want to replicate the effort in QGIS.


texts: Possibility to use spatial filters for layers in QGIS

texts:
I like the power of QGIS regarding setting attribute based filtes.

But is it also possible to filter contents of layers based on their location and spatial relations to other layers within qgis projects?

Here is an example: I only want features of one shapefile within one qgis project to be displayed (on map canvas and within the attribute table) when they touch, intersect or are disjoint from features of another shapefile within the same QGIS project.

I am not talking abaut spatial selects or exports but about instant and dynamic solutions.

I know this is supposed to be possible via using the pluging reffunctions, but it is extremely laggy and produces inconsistent results. I would love to use spatial filters outside symbology panel (rule based) as global filters


texts:
I am using osm2po to populate a pgrouting-compatible table in a postgis enabled database.
I changed the table definition a little bit (unlogged table for faster inserts, added indexes on osm_id for faster lookups, clustered on way geometry for faster seeks, ...).
I now need to update this table with fresh data. Can I somehow tell osm2po to not include DDL statements to create the table ?
Also, can the INSERT statements be changed to COPY for faster inserts ?

My current solution is to tail -n +13 osm\osm_2po_4pgr.sql to strip the first rows. It's OK, but not the best solution.
As for the COPY statements, I have no solution yet (replacing in the file would not be worth the computing time I believe).


texts: OpenStreetMap display as simple as Google Maps?

texts:
I'm glad OpenStreetMap exists at all, but I find their display has way too much details.

Here's an exemple of the same location:

Is there a way to view OSM contents with the same level of details as Google Maps, especially on smartphones?

Thank you.


texts: Calculate lat/long as percentage distance between two points (MSSQL)

texts:
I've loaded the US road shapes into MSSQL. I've determined the query to locate the segment of road for a specific address. The segment contains a "linestring" from a series of points.

I'm trying to estimate the lat/long of the actual house number.

My approach is to compare the house number to the range of numbers (so HN 13 is 30% between a road segment for house numbers 10 through 20), then determine the point of the same length along the path (so 30% down the segment, based on the earlier example)

So far i've torn the segment down into each line, and created a running total of the distance (since the house number may not be on the first or last segment)... at some point I need to identify a point along the final line... unsure how to go about.

Every calculation I see is for calculating midpoint, but given that I don't likely want a point exactly 50% between the start/end, the calculation needs to be adjusted... just not sure how.

In the unlikely event that you know of a built-in method to do this, awesome... based on my looking, I didn't find anything, so i'm not expecting.

edit: for the bounty, here's my simple test case... quick google search found the example of traveling between CYYZ (toronto) and VHHH (hong kong) to demonstrate that the process of projecting the sphere onto a 2d plane/graph doesn't always make a direct path appear straight... such an example should cover the curvature requirements (so not doing a simple algebraic lat long percentage)

DECLARE @CYYZ geography = geography::STGeomFromText('POINT(-79.6270084 43.6777176)', 4326)
DECLARE @VHHH geography = geography::STGeomFromText('POINT(113.9162921 22.3080519)', 4326)
DECLARE @path geography = geography::STLineFromText('LINESTRING(-79.6270084 43.6777176, 113.9162921 22.3080519)', 4326);

SELECT * FROM (
          SELECT 'PATH' , @path
UNION ALL SELECT 'START', @CYYZ.STBuffer(500000) -- buffer converts point to circle with
UNION ALL SELECT 'END'  , @VHHH.STBuffer(500000) -- radius, just for visibility
) A(Term, Geom)

the example is for MSSQL, since that's my platform, but the coords should work either way (note that i may have inadvertently swapped the lat/long for both points)


texts: Conflicts between localhost and Local IP address - OpenLayers

texts:
I'm having some troubles trying to visualize a tile-based map using OpenLayers.

My whole application (html, css, js and maptiles) is deployed on VirtualBox. But, when I try to access from my local machine to the map using a web browser, the map doesn't appear.

The JavaScript file contains the following lines:

var tiles = new ol.layer.Tile({
    title: 'Natural Earth',
    type: 'base',
    source: new ol.source.XYZ({
        url: 'http://localhost:8080/geoserver/www/ol_map3/{z}/{x}/{y}.png',
        attributions: 'Made with Natural Earth and TilesXYZ'
    })
})

var map = new ol.Map({
    target: 'map',
    layers: [tiles],
    view: new ol.View({
        center: ol.proj.fromLonLat([-4.72, 41.66]),
        zoom: 5
    }),
    renderer: 'canvas'
});

However, if I change the url of the tiles' source, using the local IP (192.168.56.101), the map appears on my local machine.
I don't understand what's happening...

Does anyone have an idea about the problem? Why can't I use 'localhost' instead of 192.168.56.101?


texts: ERROR 999999 during clip, Coordinates or Measures are out of bounds

texts:
The Code it occurs on each of the following interchangable lines that I have tried:

arcpy.Clip_analysis(dvof_input, feature_name, rwy_dvof_output_prefix + row['LOW_IDENT'] + rwy_dvof_output_suffix, "")
arcpy.Clip_analysis(dvof, feature_name, rwy_dvof_output_prefix + row['LOW_IDENT'] + rwy_dvof_output_suffix, "")

Note: The only difference between these two input features is one is that dvof is a clip of dvof_input occurring earlier in the script.

The Error:

ERROR 999999: Error executing function. The coordinates or measures are out of bounds. Failed to execute (Clip).

Will not work in Script, works fine when I use the Clip Tool. I even copy/pasted the exact Python Snippet from the Clip tool and put it into the Arcpy Script and the error persists.

It's even more confusing, because the same input feature clips to a different feature successfully earlier in the script (although it usually takes ~8 minutes):

arcpy.Clip_analysis(dvof_input, buffer, dvof, "")

I have another dataset of Points, similar to this dvof_input, called dof_input, which works just fine for both parts of this script:

arcpy.Clip_analysis(dof_input, buffer, dof, "")
...
arcpy.Clip_analysis(dof, feature_name, rwy_dof_output_prefix + row['LOW_IDENT'] + rwy_dof_output_suffix, "")

Worth noting, the "feature_name" polygon was created in ArcPy, and was projecting incorrectly until I added the SR to the creation of it:

## splay is a list of arcpy.Point objects ##
arcpy.CreateFeatureclass_management(arcpy.env.workspace, feature_name, "POLYGON",
                                        spatial_reference=arcpy.env.outputCoordinateSystem)
cursor = arcpy.da.InsertCursor(arcpy.env.workspace + '\\' + feature_name, ['SHAPE@'])
array = arcpy.Array(splay)
splay = arcpy.Polygon(array, arcpy.env.outputCoordinateSystem)
cursor.insertRow([splay])

Any other information I can add to diagnose this one? I tried to be comprehensive here so someone could spot a glaring issue.


texts: View ArcGIS 10 metadata in SharePoint 2010

texts:
Does anyone know how I can view the "FGDC CSDGM"-compliant metadata in SharePoint 2010?
I know I can copy the metadata output since it's an html file but that's not what I want. I want SharePoint 2010 to be able to view the metadata, or GENERATE it itself, if possible, using the xsl/xslt stylesheets in ArcGIS 10. In other words, I want to "plug in" SharePoint with ArcGIS/ArcCatalog to see 'real-time' metadata.

So this leads to my other question, does anyone know the process by which ArcCatalog 10 uses different stylesheets (there are so many) to ultimately create the final output (dhtml) which we see in ArcCatalog?

Thanks!


texts: How to get 1980 census block shapefile from 1992 Tiger/Line files

texts:
I am trying to generate a shapefile of 1980 census blocks using the 1992 Tiger/Line files. I am able to open the files directly using QGIS, and I can create shapefiles using the GDAL TIGER driver thanks to a previous answer on this forum.  However, those methods do not seem to extract the 1980 block information.  Perhaps the GDAL Tiger driver doesn't pull that info?  Is there a way to adapt these methods, or are there any other methods, for getting a 1980 census block shapefile?
Specifically, I need a 1980 block shapefile to overlay with a more recent census block shapefile to examine changed boundaries.


texts:
I am trying to save all the layers in Layers panel as MapInfo TAB file Then closed them all. But, Save as part doesn't work. There is no error but it doesn't save it. When I change the format to another format like "ESRI Shapefile" the code works well.

def packTables():

    from qgis.core import *
    from PyQt4.QtCore import *
    import os.path
    gda94CRS=QgsCoordinateReferenceSystem("EPSG:38356")
    for layer in iface.mapCanvas().layers():
        layerName=layer.name()
        layerCRS=layer.crs()
        layerSource=layer.source()
        destination="D:/1"  + "/"+layerName
        QgsVectorFileWriter.writeAsVectorFormat(layer, destination ,'System', gda94CRS,"Mapinfo TAB")
        QgsMapLayerRegistry.instance().removeMapLayer(layer)

packTables()


texts: Fix distortions when reprojecting cshapes world maps?

texts:
I am using the cshapes package in R as a source of historically accurate country borders. While the default WGS 84 maps plot fine, I would like to reproject the maps. This introduces some distortions that I'm guessing are related to underlying issues in the map data:

How can I correctly reproject and plot a cshapes map in R?
If the distortion issue really is related to underlying map data issue, how can I fix these?

Here is a default WGS 84 map, which looks fine:

But when attempting to plot a Robinson projected map, distortion near 180 longitude:

The associated shapefile is here.

I created these maps in R using the following code:

library("cshapes")
library("ggplot2")
library("rgdal")

# Get map data and plot with default CRS (WGS 84)
wmap <- cshp(date=as.Date("2012-06-30"))
wmap_df <- fortify(wmap)

ggplot(wmap_df, aes(long,lat, group=group)) +
  geom_polygon() +
  labs(title="World map (longlat)") +
  coord_equal()

ggsave("~/Desktop/map1.png", height=4, width=7)

# Reproject to Robinson
wmap_robin <- spTransform(wmap, CRS("+proj=robin"))
wmap_df_robin <- fortify(wmap_robin)
ggplot(wmap_df_robin, aes(long,lat, group=group)) +
  geom_polygon() +
  labs(title="World map (robinson)") +
  coord_equal()

ggsave("~/Desktop/map2.png", height=4, width=7)

# Export original map data
writeOGR(wmap, ".", "wmap_wgs84", driver="ESRI Shapefile")

Additional info:

Duplicate of original question on SO, someone over there suggested I try to get help here.
I know there are many other sources for country borders, but I'm tied to cshapes for historically accurate maps.
I'm using GDAL 1.11.3, installed via homebrew, and Proj4 4.9.2. I built rgdal from source against these.
Part of my problem is that I don't know what terms to even search for for this issue, any help would be great.

There is a similar question here, which indicates that the problem is that Russia and a Pacific island exceed 180°E.  However, the answers provided are either implemented in QGIS, not R, or are based on removing vertices, which is problematic for me because I want to eventually use the maps to make thematic maps (i.e. countries with color fill).


texts:
I haven't found any good explanation yet about how to display a map in MapServer (localhost).

I have ms4w with p.mapper in my pc. I have created a base .map file, shape files and html file.

Please tell me the exact steps in which folder I put these files to be displayed in p.mapper?


texts: Can group layers be published to ArcGIS Server as a single object?

texts:
ArcGIS Desktop/ArcGIS Server 10.7.1:
I have 200 layers that I generate via ArcPy on a weekly basis (from an external source).

The layers are grouped by group layers (roughly 10 layers per group).

My team plans to serve up the layers to an external system via ArcGIS Server.
However, I don't think it's wise to manually publish all 200 layers each week. There would be too much room for error and it would take too much time.
As someone who is new to ArcGIS Server/web GIS, I'm wondering,  is it possible to publish group layers as a single object, rather than as individual layers?

That would allow my team to publish 20 objects/layers each week, instead of 200 objects/layers. I think that would be much easier to manage.


texts: Moving chart to top of pop-up window in ArcGIS Online?

texts:
Is it possible to move a chart/graph to the top of the pop-up window in ArcGIS Online via its GUI?

I have numerous fields to display and it works best if the chart is displayed first.

Just the standard chart that the web map provides now would suffice.

Having to modify HTML code seems a little much for functionality that seems fairly basic.


texts: ogr2ogr does NOT save "styles", say "color", information; how do I get around it?

texts:
I'm using ogr2ogr to convert KML to GeoJSON and I'd like to extract style information along with the data.

The bash script I am using is:

filekml="/tmp/poly.kml"
filegeojson="/tmp/poly.geojson"

ogr2ogrProg="/usr/local/fwtools/FWTools-linux-i686-3.1.0/bin_safe/ogr2ogr"

# randomNum=14089918888
randomNum=`date '+%s'`

# Grab the National Weather Service Polygon (KML) data file, using SH script in Linux:
getFromURL="http://www.srh.noaa.gov/data/radar/poly.kml?${randomNum}"
wget ${getFromURL}
echo y | mv poly.kml\?${randomNum}  ${filekml}

#-- use OGR program above to convert the downloaded KML file to geoJSON format
rm -f ${filegeojson}
$ogr2ogrProg -f "GeoJSON"  ${filegeojson}  ${filekml} -nln sr_polygon

# ---- examine the KML file and its' sister geoJSON file
vi /tmp/poly.kml /tmp/poly.geojson

You will notice KML contains lines like these:

==> color  ffa9a9a9
==> scale  0.5
==> href   http://www.srh.noaa.gov/ridge/graphics/FFW.gif

BUT--- geo JSON file has absolutely NO information on STYLE (that was embedded in KML).


texts: Change GeoServer's Home Directory with Jetty

texts:
The title says it all. I would like to have GeoServer running on another path other than http://root_path:8080/geoserver/web.
I would like it to be http://root_path:8080/sub_path/geoserver/web.

The idea is having clustered GeoServers.

One running on

http://root_path:8080/sub_path1/geoserver/web

and another on

http://root_path:8080/sub_path2/geoserver/web

I am using the embed Jetty web server.
If this is not possible with Jetty, is it possible with Tomcat?


texts: Find all points within 5m with same name on large dataset

texts:
I have a large dataset (500,000 rows) in cartoddb which contains data like this.(abbreviated for easier typing)

cartodb_id | count | lat               | long                | name |  the_geom
1          |1      | 52.4118634358102  | -4.088212251663207  | aber |  theGeom
2          |1      | 52.48842511440679 | -3.4982872009277344 | card |  theGeom
3          |3      | 52.41192887963411 | -4.088115692138672  | aber |  theGeom

What I would like to do is run a query which goes through the data and for each row, finds all the points within 5 meters AND that contain the same name.

This is from a transcription project and a lot of people have clicked around the same place but not on the same place. The count highlights how many times a place has been clicked/inputted.

I would then like to create a new dataset which contains a single location and a more accurate count, as you can see from above, there are two 'aber', which are close together. For this I will just take the first lat/long and discard all the others. So I end up with a table that has a lat long as the first one but with a count of 4.

Never having played with PostGis before I seem to need to use the <> to find the closest points and also st_Distance with ::Geography. I seem to be able to find the distance from a single point using the below

SELECT *,
 ST_Distance(
    the_geom::geography,
    CDB_LatLng(53.19100611356759,-3.8425111770629883)::geography
    ) / 1000 AS dist

 FROM gb1900
 Order by dist ASC

Im just not sure how to turn this into what I need.


texts: Paste is greyed out - ArcGIS 10.4.- Can't edit in ArcMap!

texts:
I'm having a silly problem in edit session within ESRI ArcMap 10.4.
I'm trying to copy from a polygon feature to a polygon feature and when I'm using copy function the paste function turns instantly grey and can't use it. I've tried deselecting the feature, shut-down and restart of the program but with no effect.

The solution is a total restart of the computer but I don't think this is right.
Is there any other solution or someone who has/had the same problem?


texts: Problem with supervised classification (using GRASS)

texts:
So I was playing with supervised classification with GRASS. The first issue that I ran into was right after the first two steps of classification (v.to.rast and i.gensig) which was:

  unable to read signature file

Now that is an issue I've been searching the last 2 days. Is there a practical answer to that?


texts: How to change label names to numbers in QGIS

texts: I'm trying to change label names to just numbers inside each polygon. Also I would like to add a legend with the number = the name of each polygon


texts:
Can Python toolboxes be used with ArcObjects GeoProcess.AddToolbox?

I'm trying to execute a tool in a Python toolbox (.pyt) from .NET.

First I call:

GP.AddToolbox(@"C:\path\to\MyToolbox.pyt")

And that call doesn't get an error. Then I build the parameter list and then I execute it like this:

result = (IGeoProcessorResult2)GP.Execute("ToolName", parameters, null);

This immediately throws the "Error HRESULT E_FAIL has been returned from a call to a COM component." error, so nothing is returned to result.

My next step is to create a dummy script and build it into a .tbx toolbox and leave the .NET code the same to see if it's something else, but it did make me wonder if this is even possible. I couldn't fine any examples of using AddToolbox with a .pyt file.


texts:
I am trying to find a tool to create a skeleton of a large vector dataset using QGIS and its bundle. It should be directly programmatically callable through Python or commandline (QGIS, GRASS, OGR etc. are all okay).

The target is a large thin polygon with holes representing street areas and I need to obtain an approximation of the road network via skeletonization.

What I found so far:

GRASS 7.0 has v.skeleton which should be good enough but does not seem to be in for integration with QGIS in the near future.
CGAL has skeletonize algorithms, of course, but the Python support is very poor.
pySkeleton does not work for me and I probably don't have time to salvage the code.

Other alternatives (openJump skeletonizer, this old ArcGIS plugin...) cannot be easily plugged into QGIS/Python.

Any other suggestions?


texts: Unable to establish SSL connection when downloading MODIS tiles in R

texts:
Using R I've successfully used NASA's MODIS to download NDVI data for two studies, but now I'm getting this error repetitively:

    ############################
SYSTEM_WGETRC = c:/progra~1/wget/etc/wgetrc
syswgetrc = C:\PROGRA~2\GnuWin32/etc/wgetrc
--2020-04-24 08:31:46--  https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.006/2019.11.17/MOD13Q1.A2019321.h09v04.006.2019337231330.hdf
Resolving e4ftl01.cr.usgs.gov... 152.61.133.130
Connecting to e4ftl01.cr.usgs.gov|152.61.133.130|:443... connected.
Unable to establish SSL connection.
Remote connection failed! Re-try: 18
Getting file from: LAADS
############################
SYSTEM_WGETRC = c:/progra~1/wget/etc/wgetrc
syswgetrc = C:\PROGRA~2\GnuWin32/etc/wgetrc
--2020-04-24 08:31:47--  https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/6/MOD13Q1/2019/321/MOD13Q1.A2019321.h09v04.006.2019337231330.hdf
Resolving ladsweb.modaps.eosdis.nasa.gov... 198.118.194.40
Connecting to ladsweb.modaps.eosdis.nasa.gov|198.118.194.40|:443... connected.
OpenSSL: error:1407742E:SSL routines:SSL23_GET_SERVER_HELLO:tlsv1 alert protocol version
Unable to establish SSL connection.
Remote connection failed! Re-try: 19
Getting file from: LPDAAC
############################

This is a sample of what I am using to download MODIS data:

getProduct("MOD13Q1")
getProduct("MYD13Q1")
print("begin download")
tfs <- MODIS::runGdal(product = "MOD13Q1", collection = "006"
               ,tileH = 9, tileV = 4
               ,begin = "2020.01.01"
               ,end = "2020.04.10"
               ,job = "temporalComposite" # prevents new directory with timestamp from being created
               #,SDSstring = "100000000010"
               ,outDirPath = "./MODIS3/MOD13Q1.006"
               ,overwrite = TRUE)

It took about an hour to acquire data for two different tiles successfully. But then it took 26 hours to acquire data for another AND when downloaded the data was corrupt. Anyone familiar with SSL networking that can explain what went wrong?


texts: How to create azimuth raster layer by pyqgis?

texts:
I want to create an azimuth raster layer base of single point shapefile, for this I write below code that create a memory point layer that each point in memory layer is center of each pixel and have a field by name of "Azimuth", in qgis when run this program it's ok, but when I want to load created raster layer "Azimuth_lyr.tif", display this message:

What is error(s) in this code?

my code:

def azimuth_raster(self):
    # Chosen Point layer with GDAL
    pnt_layer = self.dlg.ui.comboBox_sourcelayer.currentIndex()
    pntId = self.dlg.ui.comboBox_sourcelayer.itemData(pnt_layer)
    pntLayer = QgsMapLayerRegistry.instance().mapLayer(pntId)
    feat = pntLayer.getFeatures().next()
    pntxy = feat.geometry().asPoint()

    def pixel2coord(x, y):
        """Returns global coordinates from pixel x, y coords"""
        xp = (g_pixelWidth * x) + g_originX + (g_pixelWidth/2)
        yp = (g_pixelHeight * y) + g_originY + (g_pixelHeight /2)
        return(xp, yp)

    # Open tif file
    ds = QgsRasterLayer("C:\Temp_NoiseOnePoint\prox_Distance.tif","prox_Distance")
    global g_pixelWidth, g_pixelHeight
    g_pixelWidth = ds.rasterUnitsPerPixelX()
    g_pixelHeight = ds.rasterUnitsPerPixelY()

    # extent of the layer
    ext = ds.extent()
    global g_originX ,g_originY
    g_originX ,g_originY = (ext.xMinimum(),ext.yMinimum())
    cols = ds.width()
    rows = ds.height()

    #Create temporary vector layer and add to map
    vl = QgsVectorLayer("Point?crs=", "temporary_points", "memory")
    vl.dataProvider().addAttributes([ QgsField("X_coord", QVariant.Int),
                                      QgsField("Y_coord", QVariant.Int),
                                      QgsField("Azimuth", QVariant.Int) ])
    fet = QgsFeature()
    for row in range(0,rows):
        for col in range(0,cols):
            rspnt = pixel2coord(col,row)
            fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(rspnt[0], rspnt[1])))
            fet.initAttributes(3)
            fet.setAttribute(0, rspnt[0] )
            fet.setAttribute(1, rspnt[1])
            az = pntxy.azimuth(QgsPoint(rspnt[0], rspnt[1]))
            if az < 0.0:
                azim = 360 + az
            else:
                azim = az
            fet.setAttribute(2, azim)
            vl.dataProvider().addFeatures( [fet] )
            vl.updateExtents()
            vl.updateFields()

    QgsMapLayerRegistry.instance().addMapLayers([vl])

    # Convert azimuth point layer to azimuth raster layer

    mem_layer_provider = vl.dataProvider()

    srs = osr.SpatialReference()
    srs.ImportFromWkt(vl.crs().toWkt())
    out_rst = gdal.GetDriverByName('GTiff').Create("C:\Temp_NoiseOnePoint\Azimuth_lyr.tif",
                                                   cols, rows, 1, gdal.GDT_Byte)
    out_rst.SetGeoTransform(g_x_min, g_pixelWidth, 0, g_y_max, 0, -g_pixelHeight)
    out_rst.SetProjection(srs.ExportToWkt())
    band_rst = out_rst.GetRasterBand(1)
    band_rst.SetNoDataValue(g_NoData_value)
    band_rst.FlushCache()

    gdal.RasterizeLayer(out_rst, [1], mem_layer_provider, options = ["ATTRIBUTE=Azimuth"])

    # close bands
    band_rst = None

    # close rasters
    out_rst = None


texts: ArcGIS Pro interpreting Viewshed results - number of Observers

texts:
This question seeks to clarify a point made in an answer to ArcGIS Viewshed without Observer Points, which states:

I'm fairly sure that the number of points [Viewshed] generates along the line depends on the DEM resolution; for example, if you've got a 100m line and a raster with 10m cells, your line will be analyzed as 10 points.

Can anyone confirm this is true in ArcGIS Pro? I have a DEM raster of 30-meter cell size, and a polyline observer feature that is quite long (a 463 kilometer boundary line to be precise.) I'm getting a maximum of 279 observers. When you work out the math that would mean ≤2% of the boundary can "see" any given cell. It seems too low.
(Specifically I'm looking to convert the number of observers into the percentage of the boundary.)


texts: How to create service url to get demographic data and then tie up the demographic data with zip code?

texts: I want to create a Map tool like https://www.eddmexperts.com/eddm-map-tool.  This site is using ArcGIS map service but I did not understand where I can get the service URL they are using to get demographic data.


texts: how to create the Layer Properties button as geoexplorer?

texts:
I was in the debug mode.And I create a Layer Properties button ,but it is only show the title.
the code for it,

 * @require plugins/LayerProperties.js
--------------------
{
  ptype: "gxp_layerproperties",
  actionTarget: ["tree.tbar"]
}
-----------------------

How can I do to get a button as Layer Properties on geoexplorer tree panel that has 3 tags,like about,display and available styles? Especially the display tag is more important to move opacity slidebar.


texts: Update carto map with new dataset without re-styling

texts:
I have been using a dataset "dataset 1" for a Carto map for some time, so it is now outdated. I have an updated dataset that contains more row data but same structure.

How do I update my map with the new dataset without having to re-style the map again? I am using Carto builder.

Here Can I change an existing Cartodb map to use a different dataset?, the user offers a solution, which is create a new layer using the updated dataset, then copy and paste the previous CSS into the new layer.

However, I would still need to re-style the pop-up info window. Furthermore, will I also have to re-create my widgets?


texts: QGIS SVG icons in print preview are showing up as white silhouettes

texts: I have created my own SVG icons in Illustrator. They show up fine in the working QGIS map, however when I create the print preview, many, but not all of the SVG icons loose layers, and most show up as a white silhouette.


texts: Progressive union of a large number of multiband rasters in PostGIS

texts:
I am deriving multiband rasters from airborne lidar imagery, and loading these into PostGIS along with associated metadata. Raster cell values are counts of certain point classes over a series of vertical strata. Each input raster is 200 x 200 cells (with slight variation) and has 160 bands. Rasters from adjacent images overlap slightly. Stitching them together with ST_Union and the 'SUM' option works perfectly when dealing with small numbers of rasters, but I've been stuck for a while on how to scale this up to the many thousands to be processed without things slowing to a crawl. I am relatively new to PostGIS, so I'm looking for advice about table design and query structure (detailed below).

The current design involves two raster tables. The first table (A) holds rasters as derived from each lidar image, with no regard to overlap. Tile sizes in this table are taken from the inputs and can vary slightly. The second table (U) holds a unioned version of the same data with a uniform 100x100 tile size. This duplication of data across two tables is motivated by the need to correct and re-import data when problems are identified in individual source LiDAR images, something that occurs commonly.

Importing each new raster involves the following steps:

Load raster into A. Tile size is taken from input and can vary between recs.
Identify any 100x100 raster tiles in U that overlap the new raster in A.
In a temp table (T), union (SUM) the new raster with the overlapping tiles from U.
Delete the overlapping tiles from table U.
Insert the unioned data from T into B.

This seems to work but slows to a snail's pace as the number of rasters grows. There is a spatial index on table U and the number of overlapping rasters cannot be more than a few by definition, so I expected the running time for each iteration of the above steps to be more or less constant. Obviously I'm missing something.

As an experiment, I imported larger batches of rasters (several hundred) into table A and then attempted to union these in one go, but after more than 12 hours this had still not completed.

The initial processing of the lidar images and subsequent import into PostGIS is done via R using custom written R packages. The steps detailed above are performed by sending SQL from R to PG. I've experimented with running the same queries via psql but performance appears identical, so it doesn't seem to be a problem with R per se, more likely a problem with my SQL code.

Here's the code:

-- Identify existing rasters that overlap the new data
create or replace view mgazone56.overlaps as
  select pcu.rid from
    mgazone56.point_counts_union as pcu, mgazone56.tmp_load as tl
    where st_intersects(pcu.rast, tl.rast);

create table if not exists mgazone56.tmp_union (rast raster);

truncate table mgazone56.tmp_union;

-- Union overlapping rasters with new data, summing overlap values
insert into mgazone56.tmp_union
  select ST_Union(r.rast, 'SUM') as rast from
    (select rast from mgazone56.point_counts_union
       where rid in (select rid from mgazone56.overlaps)
     union all
     select rast from mgazone56.tmp_load) as r;

-- Delete overlapping rasters from main table
delete from mgazone56.point_counts_union
  where rid in (select rid from mgazone56.overlaps);

-- Insert updated data
select DropRasterConstraints('mgazone56'::name, 'point_counts_union'::name, 'rast'::name);

insert into mgazone56.point_counts_union (rast)
  select ST_Tile(rast, 100, 100) as raster
  from mgazone56.tmp_union;

select AddRasterConstraints('mgazone56'::name, 'point_counts_union'::name, 'rast'::name);

-- Delete records from the temporary tables
truncate table mgazone56.tmp_load;
truncate table mgazone56.tmp_union;


texts: Where can I make a modification or suggestion for Google Geocode API response, specifically the formatted_address field?

texts:
I was wondering where I can make a suggestion to modify a Google Maps geocode API response. Specifically the formatted_address field

Here is why : Usually airport type locations always have in their formatted_address field the airport code e.g: (JFK). This is the case for almost every airport I tried to geocode.

But this one Aéroport Marrakech-Ménara is missing the airport code

Where do you think I can suggest this modification so that the formatted_address follows the standard?

Thanks in advance.


texts: Moran's I test on regression residuals with spatial autocorrelation

texts:
Suppose we want a regression model $y_i = f(x_i^t) + e_i$, where $y_i$ and $x_i$ are collected spatially, and we believe there would be spatial autocorrelation (SAC). I read some popular methods like spatial-lag or spatial-error, and the Moran's I test.

My question is, if we believe that the SAC in $y_i$ is all caused by SAC $x_i$, that is, we believe $e_i$ are independent with each other. Then is there any problem to simply run a normal regression model? It looks to me that this way does not violate any assumptions. I can definitely run a test on the residuals. However several application papers do things in such a way: they first run Moran's I on $y_i$, if significant, then spatial-lag. O.w., run Moran's I on residuals, if significant, then spatial-error. If looks weird to me.


texts: What is the more accepted practice for defining tolerance values?

texts:
The Context
Due to restrictions, I have an application that requires line features to be selected if they are within a user defined tolerance from a points data set.  I can not currently use any third party API to perform so have written an efficient algorithm that completes this task effectively.  The application will be used by people that may or may not have a sufficient grasp on spatial concepts such as the defined unit of measure for their data set, or parameters requiring 'acceptable' values in map units.

My Befuddlement

The use of the Tolerance parameter and whether I should keep the tolerance in a commonly known unit such as Centimetres, Meters or Kilometres and perform the required calculations to massage the input into an appropriate value; or
Assume the user knows the data, or has the ability to ascertain the data sets map units and enter an intelligent, useful parameter value.

The Question:
What is the more accepted practice?

The Implications
If Befuddlement 1 is the accepted practice or considered more user friendly, can you proved the some direction towards ascertaining the Datum and other relevant coordinate information to perform the required function.

If Befuddlement 2 is more generally accepted, then case solved.  Though it does make the solution easier and more straight forward, I feel forcing the user to go through a discovery process may lessen the user experience of the software.

Any thoughts on the matter would be greatly appreciated.


texts: How to add multiple labels and points in OpenLayers 3 using WMS layer

texts: I have India country WMS layer, and I want to add label for each state and point or small circle for each state. I created a layer from GeoServer using a PostgreSQL database. Using the database how can I show the labels and points for each state with OpenLayers 3?


texts: error of data conversion value getting while importing shapefile in shp2pgsql-gui on ubuntu

texts:
Connecting:  host=localhost port=5432 user=postgres dbname=tremblement_de_terre
Connection succeeded.
Connecting:  host=localhost port=5432 user=postgres dbname=tremblement_de_terre

==============================
Importing with configuration: m_cap point, public, geom, /home/jean/Documents/M_CAP Point.shp, mode=c, dump=1, simple=1, geography=0, index=1, shape=1, srid=0
Shapefile type: MultiPoint
PostGIS type: POINT[2]
Unable to convert data value to UTF-8 (iconv reports "Chaîne multi-octets ou étendue de caractères invalide ou incomplète"). Current encoding is "UTF-8". Try "LATIN1" (Western European), or one of the values described at http://www.postgresql.org/docs/current/static/multibyte.html.

Shapefile import failed.

After trying to know the encode setting with the ubuntu shell command  file -f .shp,the result is something like

file -f M_CAP\ Point.shp
:                ERROR: cannot open ' (No such file or directory)
:                ERROR: cannot open' (No such file or directory)
:                ERROR: cannot open ' (No such file or directory)
��OA�\���?A�: ERROR: cannot open��OA�\���?A�' (No such file or directory)
��OA:         ERROR: cannot open ��OA' (No such file or directory)
��OA:          ERROR: cannot open��OA' (No such file or directory)
�5A�X=�:         ERROR: cannot open �5A�X=�' (No such file or directory)
�5A�X=�:         ERROR: cannot open�5A�X=�' (No such file or directory)
�5A�X=�:         ERROR: cannot open �5A�X=�' (No such file or directory)
��XCAp���_A~�: ERROR: cannot open��XCAp���_A~�' (No such file or directory)

When I also try to open this file with QGis as a Vector layer, using UTF-8 encoding, It opens clearly. What happens now when trying to import it to postgres using shp2pgsql-gui?

Please help me, the original file have been created on windows, and I am using Ubuntu 12.04 Lts. Thanks a lot


texts: What is the meaning of DN and DRVAL2 categorie?

texts:
I try to find the meaning of the categories DN and DRVAL2 when I used them in the command:

lrenderer = new QgsCategorizedSymbolRendererV2("DRVAL2" , Categrories );

Is there any tutorial or documentation about all of this kind of parameter ?


texts: In a webgis using leaflet, how to export maps as pdf?

texts:
Like the title says, I'm building a webgis using Leaflet. I want to allow users to export maps from the website and download them, but the current downloaded format is shapefile. I wonder if there's a configuration that enables export format as pdf.

As it may be related, the maps only show the regions in which data are visualized and hide the rest of Leaflet base map. Can you tell if this configuration will be preserved in pdf exported outputs?


texts: GDAL RasterizeLayer does not Burn All Polygons to Raster?

texts:
I am trying to burn a shapefile to a raster using GDAL's RasterizeLayer. I pre-create an area of interest raster from a different shapefile, given a specific pixel size. This AOI then serves as a base for all following rasterizations (same number of collumns and rows, same projection and geotransform).

The issue occurs, however, when I go to burn the shapes to their own raster, based on the same pixel size and projections. The link below (don't have enough rep to post the image), shows the original shapefile in tan, and the dark pink where RasterizeLayer has burned data. The light pink is nodata values for the dark pink raster data. The grey is the AOI based on which the shapefile burn was completed.

Given the extents of the shapefile polygons, I would expect to see burn values in the bottom two corners, as well as the two pixels underneath the data that does show. Obviously, however, this is not the case.

As follows is the code that I used to generate these. All of the shapes were created using QGIS, and were all created in the same projection.  (It should be noted that the gridding in the picture shown was just to give an idea of the pixel size I was using.)

from osgeo import ogr
from osgeo import gdal

aoi_uri = 'AOI_Raster.tif'
aoi_raster = gdal.Open(aoi_uri)

def new_raster_from_base(base, outputURI, format, nodata, datatype):

    cols = base.RasterXSize
    rows = base.RasterYSize
    projection = base.GetProjection()
    geotransform = base.GetGeoTransform()
    bands = base.RasterCount

    driver = gdal.GetDriverByName(format)

    new_raster = driver.Create(str(outputURI), cols, rows, bands, datatype)
    new_raster.SetProjection(projection)
    new_raster.SetGeoTransform(geotransform)

    for i in range(bands):
        new_raster.GetRasterBand(i + 1).SetNoDataValue(nodata)
        new_raster.GetRasterBand(i + 1).Fill(nodata)

    return new_raster

shape_uri = 'activity_3.shp'
shape_datasource = ogr.Open(shape_uri)
shape_layer = shape_datasource.GetLayer()

raster_out = 'new_raster.tif'

raster_dataset = new_raster_from_base(aoi_raster, raster_out, 'GTiff',
                                -1, gdal.GDT_Int32)
band = raster_dataset.GetRasterBand(1)
nodata = band.GetNoDataValue()

band.Fill(nodata)

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, burn_values=[1])

Is this a bug in GDAL, or is RasterizeLayer burning data based on something other than just the presence or lack of a polygon within a specified pixel area?

The files that I was using can be found here.


texts: Creating boundaries around lines based on attribute

texts:
I am trying to create polygons around lines within another polygon based on a lines attribute

Basically from this I want three boundaries within this boundary for each section of line. I have been scouring for solutions but have yet to find one.

Unfortunately I only have an ArcGIS standard license.


texts:
I have to create vegetation and sensitivity maps for sites. I have a site boundary shapefile and map from raster aerial imagery behind the site boundary vector.

When I used to use Arc I used the cut polygon tool. This allowed me to create a one shapefile vegmap with different polygons that I could label according to what they are (eg: forest, grassland, disturbed areas). The cool thing about this is that it eliminated any of those slivers of areas not mapped or intersecting polygons). This is important as I need to use the map to calculate areas of each of the different vegtypes.

I cannot seem to find a similar function in QGIS. I have fiddled with snapping and it really isn't good enough. It will take me hours and hours to individually snap the vertices for a small map, let alone a huge and intricate one!

Does anyone have an answer for me to do this in QGIS?

I have not got ArcGIS Desktop and cannot do it there.


texts: Why does bindPopup not work inside my transparent Leaflet polygons?

texts:
I am trying to display a popup from the attribute value of a geojson layer by calling the geojson file. However, the popup displayed on the edges of the boundaries of the polygons features only.
Here is the code:
var region = L.geoJson(region_data, {
    style: {
        fillColor: 'none',
        weight: 2,
        opacity: 1,
        color: 'black'
    },
    onEachFeature: function (feature, layer) {
        if (feature.properties && feature.properties.name) {
            layer.bindPopup(feature.properties.name);
        }

    }
})


texts: Seeding Tilecache with gdal2tiles output

texts:
Very simple question. I've executed gdal2tiles against a GeoTIFF coverage, composed by several TIFF files. I want to serve the resulting tiles using TileCache.

In the TileCache documentation, I see how to create new layers, using as sources JPG images, WMS layers, MapServer layers, etc (examples here). But the point is I want to directly provide the tiles generated by gdal2tiles as source. So, I want to access my tiles with something like:

http://example.com/yourdir/tilecache.cgi?LAYERS=basic&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG:4326&BBOX=-180,-90,0,90 &WIDTH=256&HEIGHT=256

Or

http://example.com/yourdir/tilecache.cgi/1.0.0/basic/0/0/0.png

Using any TMS client. Now I'm using the Google Maps viewer generated by gdal2tiles, based on GTileLayer

Any tips on TileLayer configuration to achieve this goal are welcome. Or any suggestion to serve gdal2tiles output (tiles following the TMS spec) in a different way. I know, gdal2tiles was built following the TMS specs, before the WMTS standard. So, there're probably better ways to serve the gdal2tiles output today.

Many thanks in advance!


texts: Open popups in a fixed position with Leaflet

texts:
I have create my map with Leaflet exactly how I want it (with my markers and externals links).

I would like to have the popup always fixes in the top right hand corner (for an example).

Is there any way to fix the popup's location ?


texts: ArcMap 10.5.1 address locator properties via arcpy

texts:
I am making address locators in arcpy but need to set 2 properties values after the locator is created. see the image below.. the property values for "match without house number" and "match with no zones" both need to be set to value "Yes"

EDIT**

this cannot be accessed via Arcpy but the .loc file will contain these values.
the issue is that the .loc file originally will not contain the fields to edit the values you have to manually edit them within catalog/arcmap then the .loc file will contain the values. the code below works but only after you enable these extended properties within arcmap/catalog which defeats the prupose of automating this task, I do not want to manually go into arcmap/catalog to set the values.

with open(r"Y:\pathtolocator.loc", 'r') as file :
  filedata = file.read()

# Replace the target string
filedata1 = filedata.replace("supportsEmptyHouseNumber = false", "supportsEmptyHouseNumber = true")

filedata2 = filedata1.replace("supportsOptionalZone = false", "supportsOptionalZone = true")

# Write the file out again
with open(r"Y:\pathtolocator.loc", 'w') as file:
  file.write(filedata1)
  file.write(filedata2)


texts: Service Area polygon appears to fill in between the roads for areas it should not reach

texts:
In ArcMap 10.5.1, I have created a service area polygon of 1000m along a road network, with a trim of 50m. In the image you will see where the road is linked around, the inside is filled in (red circles). However there should be a hole in this area as it is outside of the trim and not reachable within the parameters.

I am using ArcGIS Network Analyst. After adding the Network Dataset I used the Network Analyst toolbar to create a new service area. I added Facilities, which was one single central point. Within Service Area layer properties I set the Impedance to 1000m, in the direction away from the facility. I chose to Generate Polygons - Detailed with a 50m trim. I also chose to generate a line.

What am I doing wrong here?


texts:
I have many netcdf files in one folder. I want to extract variables from each netcdf file and save them as another .nc file.

The process is like this:

Make NetCDF Raster layer (one .nc)
Define name, dimension, 2 variables
Make NetCDF from Raster
Define name, dimension, 2 variables saved in 2 different .nc files

This is okay with only few files. But I have many files, so how do I do this in model builder?

 I am using ArcMap 10.2.


texts: Convert Mapinfo files into KML files, to check if coordinates are inside a polygon

texts:
Sorry for my absolute ignorance, I have 1 day of experience in GIS systems (I'm just a Java developer).
Our client did provide us a compressed file with huge .TAB, .MAP, .ID, and .MAP files, and wants a little app to check if certain coordinates are inside the "coverage" polygons (defined on the aforementioned maps).
EDIT:
I've found that the files that I have are MapInfo files, and they can be converted into Google KML files, which can be later processed to see if certain coordinates are inside. However, the online tools I found to do this conversion are size limited (ex: 5 MB and my MAPINFO file is 245 MB)
How could I make this conversion?


texts: Limit 1 route to leaflet L.Routing.control

texts:
I'm using the leaflet-routing library https://github.com/perliedman/leaflet-routing-machine/releases/tag/v3.1.1

How can I limit the amount of routes the L.Routing.control returns.

I'm using this code:

fDirectionService = new L.Routing.control({
    waypoints: [
        new L.latLng(-27.1216666667, -52.6213888889),
        new L.latLng(-27.5963, -48.5491)
    ]
}).addTo(map).on('routesfound', function(e) {
    console.log(e.routes); // e.routes have length 2
});

After the 'routesfound' event it returns 2 routes (fDirectionService._routes.length == 2)

Can I request only one? Or at least hide all except the first?

Edit: There is an option showAlternatives: http://www.liedman.net/leaflet-routing-machine/api/ but even that doesn't work.

Edit 2: Days and nobody answered, that's odd. Anyway, I've tested with Mapbox e with the same coordinates it returns only one route, could it be that it must be configured on the server side only?


texts: QGIS: On the fly projection from WGS84, 19 N, to Web Mercator. Do I need to define a new CRS with transformation parameters?

texts:
Workin in QGIS, I have a shapefile called "protected areas" of Dominican Republic, that was digitized from topomaps 50k. Shapefile is projected to WGS84, 19N zone. My problem is that when I load openlayers, the limits of almost all protected areas, don't get aligned with openlayers. In the following image there is an example:

The limit of the protected area (in red) is the coastline; but the real coastline plotted in Openlayers (google satellite) doesn't align with my shapefile (this area has rocky coastline, that haven't changed since topo was made). There is a ca. 45m shift in the X axis, and ca. 50m in the Y axis. In the following image:

the limit of the protected area is prefectly aligned with topomap (coastline, road, etc.). As I know, OpenLayers is projected using "Web Mercator". So, do I need to define a new CRS with transformation parameters between WGS84 19N and Web Mercator? I believed they use the same ellipsoidal reference, so I shouldn't need to define a transformation.

Any help would be appreciated.

Thanks.

Jose


texts: Calculating distances between points for every available point in the dataset

texts:
The data of interest: a point vector and a line vector. Both shapefiles, with the line vector acting as the road/motorway. All of the points are along the road so close by.

I want to calculate the distance from point to point but the points are not necessarily in order in terms of their spatial representation.

What I mean by this is that point 64 is one the right side of the road, while point 65 on the left side of the road (I have attributes for that, road, side). E.g. id 65 might not be the next point after 64 so as in a straight line.

So far, I have tried to use distance matrix but that doesn't take into account the motorway direction (line vector).

Any ideas what to do next?


texts: We are using mapquest tiles(ol.source.MapQuest) in our application. But, mapquest team is discontinued to access their tiles from 11th July, is there any open source alternative to use instead of mapquest.


texts: How do I change the ArcGIS Server Web Adaptor name?

texts:
If I decide I want to change the Web Adaptor name, is it as simple as going to the Web Adapter tab in the Server Manager and clicking unregister?

And then would I just rerun the Web Adaptor install and choose a different name?

Or what is the proper way of changing Web Adaptor names?


texts: Openlayers 3 getting 401 error from geoserver with secured layers

texts:
I have an OL3 webmap that is attempting to access services from Geoserver that is running on a server on our internal network that uses LDAP logins using windows credentials.  I found this link (Accessing secure Geoserver layer with username password) and this link (https://stackoverflow.com/questions/10950747/authenticate-in-geoserver-with-asp-net-and-iis/13123640#13123640) and tried both but it is not resolving the 401 error issue.

This is the code I have from the second link:

//login function
function login (options) {
// url del servlet del geoserver
var url = options.server + "/geoserver/j_spring_security_check";
// parametros para el login
params = "username=" + options["user"] + "&password="
            + options["password"];

var contentType = "application/x-www-form-urlencoded";
//se inicializa la petición ajax
var ajax = $.ajax({
    data : params,
    type : "POST",
    contentType : contentType,
    url : url
});
// se ejecuta cuando la peticion finaliza
ajax.done(function() {

    if ($.cookie("JSESSIONID") != null && options && options.success) {
        options.success();
    }
});
// si ocurrio un error al realizar la peticion
ajax.fail(function(data) {
    if (options && options.failure) {
        options.failure(data);
    }
});
// se ejecuta siempre al final de la petición, sin importar que esta
// haya fallado
ajax.always(function() {
    if (options && options.always) {
        options.always();
    }
});
};
var un = prompt("enter your username","");
var pw = prompt("enter your password","");
login({
    user: un, //geoserver user
    password: pw,
    server : "http://(host server here):8080", //geoserver host
    success : function(){
        alert("Login OK!");
    },
    failure : function(){
        alert("Login fail!");
    }
});

Also, if I log into the geoserver admin webpage in a separate tab, I can then open the webmapping application with no issue.  I would like to ideally remove this step by putting in the geoserver login into the actual webpage (and preferably not store the login credentials anywhere for security purposes).


texts: How to limit WMS to a certain zoom level in QGIS?

texts:
I am trying to export a map to an image at 600DPI and 1:2,000,000 scale but have the background WMS map stay at a certain zoom level as the text can get really small compared to the other layers I have.

It'll keep the zoom level I want and that is being displayed if I go to export the map as image and just keep the resolution setting as 96DPI but the other layers are 600DPI and I want to keep the quality for printing.

The background layer is WMS

How can I set a certain zoom level so whenever I try to export with whatever settings it will only use that zoom level?

The WMS server is: http://gaservices.ga.gov.au/site_7/services/NationalMap_Colour_Topographic_Base_World_WM/MapServer/WMSServer

and I'm trying to have the zoom level set to look like this https://i.imgur.com/JthE9Fu.png which has nice big text and limited roads etc


texts: Error in merging two raster files (precipitation) in QGIS

texts: The newly computed geotiff shows the extent of both raster layers, but only the data of one. If I use identify features tool it gives me only the values of the visible part not of the other part. I use Float32. What am I doing wrong?


texts: Improve performance on a st_dwithin query (in PostGIS)

texts:
At the momenent I am working on a query, as described in an earlier question. I have the two tables

testme (tracked GPS-profile with Point Geometry) and
roads (geometry roads shapefile)

Besides the distance of each tracked point (row) in table pt in want to find the closest point using a combination of ST_Distance and ST_DWithin.

DROP TABLE IF EXISTS raw_2015_processed;
EXPLAIN ANALYZE
CREATE TABLE raw_2015_processed AS
SELECT pt.id,
    pt."DeliveryID",
    pt."VehicleID",
    pt."TrackID",
    pt."Longitude",
    pt."Latitude",
    pt."Altitude",
    pt."Heading",
    pt."Speed",
    pt."Satelites",
    pt."HDOP",
    pt."VDOP",
    pt."Xfcd",
    pt.ts,
    pt.received,
    pt.the_geom,
(SELECT ST_ClosestPoint(line.geom,pt.the_geom) AS closest_geom
    FROM roads AS line
    WHERE ST_DWithin(line.geom,pt.the_geom, 0.5) LIMIT 1),
(SELECT ST_Distance(line.geom,pt.the_geom) AS distance
    FROM roads AS line
    ORDER BY pt.the_geom <#> line.geom LIMIT 1),
(SELECT ST_AsText(ST_ClosestPoint(line.geom, pt.the_geom)) AS closest_coordinates
    FROM roads AS line
    ORDER BY pt.the_geom <#> line.geom LIMIT 1)
FROM raw_2015 AS pt
ORDER by pt.id;`

Working with a reduced file 'Testme' of 144 rows EXPLAIN ANALYZE returns following Query Plan:

"Sort  (cost=84207.04..84207.40 rows=144 width=131) (actual time=11797.677..11797.685 rows=144 loops=1)"
"  Sort Key: pt.id"
"  Sort Method: quicksort  Memory: 54kB"
"  ->  Seq Scan on testme pt  (cost=0.00..84201.87 rows=144 width=131) (actual time=82.458..11797.268 rows=144 loops=1)"
"        SubPlan 1"
"          ->  Limit  (cost=9.46..578.16 rows=1 width=155) (actual time=81.508..81.508 rows=1 loops=144)"
"                ->  Bitmap Heap Scan on roads line  (cost=9.46..578.16 rows=1 width=155) (actual time=81.388..81.388 rows=1 loops=144)"
"                      Recheck Cond: (geom && st_expand(pt.the_geom, 0.5::double precision))"
"                      Rows Removed by Index Recheck: 0"
"                      Filter: ((pt.the_geom && st_expand(geom, 0.5::double precision)) AND _st_dwithin(geom, pt.the_geom, 0.5::double precision))"
"                      ->  Bitmap Index Scan on geom_index_roads  (cost=0.00..9.46 rows=139 width=0) (actual time=79.722..79.722 rows=450402 loops=144)"
"                            Index Cond: (geom && st_expand(pt.the_geom, 0.5::double precision))"
"        SubPlan 2"

However, when I run the query for a larger dataset with more data points (< 5 Million points) , it gets very slow (i.e. several hours / days). Do you guys see a way to increase the speed of the query? Is there an an alternative option to st_dwithin or a different query structure which can proof to be helpful?


texts: Splitting lines at vertices in ArcGIS Desktop with Basic level license?

texts:
Is there a way to work around splitting lines at vertices in ArcGIS Desktop without an Advanced level license?

I don't have the advanced license and I have a model that uses this tool.


texts: How can i resolve problem with bounds for export?

texts:
How could I set up bounds in TileMill if i don't want to upload the country around? I simply want to upload my created island in TileMill to Mapbox.
Is there a way to do it? I am using OSM Bright for the design of my map.


texts: Find spatial and temporal variability from a GeoPandas dataframe of polygons

texts:
I have a geopandas dataframe of polygons (all rectangles ~100mx120m) and each of these polygons has an attribute for time and % ice coverage. This is all derived from imagery I collected in the Antarctic using a fixed wing drone and I have already turned the imagery into a calculation of ice coverage in each image. I now want to look at the spatial and temporal variability of ice coverage in these areas.

Ideally I would get some quantitate values telling me at what spatial distance ice coverage is no longer correlated (decorrelation scale) and after what length of time the ice coverage stops being correlated. I have thought about using variograms, but I'm not clear on how to do it with data of this form (polygons randomly spatially distributed with some spatial overlap on different surveys).

Here is what a subset of my survey looks like generally:

Any thoughts for how to do this? pysal seems promising but I'm struggling a bit with the documentation and I'm not sure if there is a better option for integrating with geopandas or just doing this without a specific package.

The files are available here and can of course be opened and inspected with the code below:

import geopandas as gpd
gdf = gpd.read_file('disco_coverage.shp')
gdf.head()


texts: Updating an Openstreetmap PostgreSQL database after the system restore

texts:
My Openstreetmap PostgreSQL database takes a long time – about 25 hours by now – for updating via Osmosis after I had made a system restore at the last weekend. Before making the system restore it took just 18-27 minutes for updating my database.

I work with Ubuntu 14.04 (Trusty) Linux operating system. The Ubuntu system and the Openstreetmap PostgreSQL database are situated on different hard drives. I ought to add that I made the system restore just for a hard drive where the Ubuntu operating system had been. I didn't touch the disk with the Openstreetmap database.

So my questions are as follows. What had I done wrongly and how may I escape such the problem in the future? Can I speed up an updating process if I will take a backup copy of my database? I made that backup copy in October.


texts: DLL As Embedded Resource ArcMap Add-In? VB.NET

texts:
I tried using StackOverflow, but I don't think anyone can answer my question in context(ArcMap Addins).

I am currently developing an add-in that needs to load a custom control DLL as a dependency. Unfortunately, ArcMap uses DLL's as their add-ins. So in turn, I must be able to either import an embedded resource to the Visual Studio Toolbox(for use with the forms), or figure out a way to make the control reference the embedded DLL as its delegate.

Let me know if you happen to have any ideas.

Thanks,

Logan


texts: Cannot get changeAttributeValues() to work in a loop

texts:
I'm trying to use QGIS to loop through a set of points, find their elevation from a DEM and use that value to update the point layer's attribute table.  The loop works fine until I try to use changeAttributeValues(), and I can change the attribute values outside of the loop, but when I try to put the two together I cause a dump.

while provider.nextFeature(feat):
    geom = feat.geometry()
    x = geom.asPoint()
    res, ident = rlayer.identify(QgsPoint(x))

    for (k,v) in ident.iteritems():
        elevation = float(v)
    attrs = { 2 : elevation }
    caps = vlayer.dataProvider().capabilities()
    if caps & QgsVectorDataProvider.ChangeAttributeValues:
        vlayer.dataProvider().changeAttributeValues({ fid : attrs })

This causes the following error in Linux:

*** glibc detected *** python2: free(): invalid pointer: 0x00000000024655b8 ***

Edit:

The problem isn't the loop, it is now I'm calling changeAttributeValues().  I've simplified the script down to just update a single point and it causes the same core dump.  Here's the updated script:

#!/usr/bin/env python2
# -*- coding: utf-8 -*-

from qgis.core import *
import qgis.utils

#initializes QGIS and points python to where it can find QGIS's stuff
QgsApplication.setPrefixPath("/usr/", True)
QgsApplication.initQgis()

vlayer = QgsVectorLayer("/gis/vector/elev_pts.shp", "elev_pts", "ogr")
if not vlayer.isValid():
    print "Layer did not load!"

vlayer.dataProvider().changeAttributeValues({ 0: 2: 432.1 })

QgsApplication.exitQgis()


texts:
I am looking for an alternative to ESRI's Simplify Polygons - Bend Simplify method (I only have a Basic License).

Is there an existing QGIS plugin which supports this specific polygon simplifying algorithm?

If not, does anyone have a suggestion for which Generalization algorithm within the Grass v.generalize might produce the most similar results?


texts: How to convert lines to polygons using NetTopologySuite?

texts:
I am trying to convert lines to polygons using NetTopologySuite(a C# port of Java Topology Suite JTS). I ran the sample console code succesfully and it generated two polylines.

I then tried to run the code using three line strings below which are noded at endpoints.

(0 0,10 0)
(10 0,10 10)
(10 10,0 0)

but this move is bringing the demons out of the code. A bunch of errors are popping up as listed below.

C:\Users\Home\Desktop\polgon\polcreate\polcreate\bin\Release>polcreate.exe
   at NetTopologySuite.Geometries.LinearRing`1.ValidateConstruction() in E:\Programming
\JTS\JTS2.11\v2.11_cleanup\NetTopologySuite\Geometries\LinearRing.cs:line 119
   at NetTopologySuite.Geometries.LinearRing`1..ctor(ICoordinateSequence`1 coordinates,
 IGeometryFactory`1 factory) in E:\Programming\JTS\JTS2.11\v2.11_cleanup\NetTopologySui
te\Geometries\LinearRing.cs:line 43
   at NetTopologySuite.Geometries.GeometryFactory`1.CreateLinearRing(ICoordinateSequenc
e`1 coordinates) in E:\Programming\JTS\JTS2.11\v2.11_cleanup\NetTopologySuite\Geometrie
s\GeometryFactory.cs:line 552
   at NetTopologySuite.Geometries.GeometryFactory`1.CreateLinearRing(IEnumerable`1 coor
dinates) in E:\Programming\JTS\JTS2.11\v2.11_cleanup\NetTopologySuite\Geometries\Geomet
ryFactory.cs:line 535

and so on....

The full code is shown below. Can anyone guide me the right way?

#define simple
using System.Text;
using System;
using System.Collections.Generic;
using GeoAPI.Geometries;
using GeoAPI.IO.WellKnownText;
using NetTopologySuite.Geometries;
using NetTopologySuite.Operation.Polygonize;

#if simple
using NUnit.Framework;
using Coord = NetTopologySuite.Coordinates.Coordinate;
using CoordSeqFac = NetTopologySuite.Coordinates.CoordinateSequenceFactory;
#else
using Coord = NetTopologySuite.Coordinates.BufferedCoordinate;
using CoordSeqFac = NetTopologySuite.Coordinates.BufferedCoordinateSequenceFactory;
#endif

namespace polcreate
{
    /// <summary>
    /// Example of using Polygonizer class to polygonize a set of fully noded linestrings.
    /// </summary>
    [TestFixture]
    public class PolygonizeExample
    {
        [STAThread]
        public static void Main(String[] args)
        {
            PolygonizeExample test = new PolygonizeExample();
            try
            {
                test.Run();
            }
            catch (Exception ex)
            {
                Console.WriteLine(ex.StackTrace);
            }
        }

        [Test]
        public void Run()
        {
            IGeometryFactory<Coord> geoFactory =
                new GeometryFactory<Coord>(new CoordSeqFac());
            WktReader<Coord> rdr
                = new WktReader<Coord>(geoFactory, null);
            List<IGeometry<Coord>> lines
                = new List<IGeometry<Coord>>();

            // old values which ran successfully
            //lines.Add(rdr.Read("LINESTRING (0 0 , 10 10)"));
            //lines.Add(rdr.Read("LINESTRING (185 221, 100 100)")); //dangling edge
            //lines.Add(rdr.Read("LINESTRING (185 221, 88 275, 180 316)"));
            //lines.Add(rdr.Read("LINESTRING (185 221, 292 281, 180 316)"));
            //lines.Add(rdr.Read("LINESTRING (189 98, 83 187, 185 221)"));
            //lines.Add(rdr.Read("LINESTRING (189 98, 325 168, 185 221)"));

            //i changed the values here..
            lines.Add(rdr.Read("LINESTRING (0 0 , 10 0)"));
            lines.Add(rdr.Read("LINESTRING (10 0, 10 10)"));
            lines.Add(rdr.Read("LINESTRING (10 10, 0 0)"));

            Polygonizer<Coord> polygonizer
                = new Polygonizer<Coord>();
            polygonizer.Add(lines);

            IList<IPolygon<Coord>> polys = polygonizer.Polygons;

            Console.WriteLine("Polygons formed (" + polys.Count + "):");
            foreach (object obj in polys)
            {
                Console.WriteLine(obj);
            }
            Console.ReadLine();
        }
    }
}


texts: How can I geocode to a shape instead of a coordinate?

texts:
Geocoding to a shape instead of a point

The project I'm currently working on is Find-A-Record. We are geocoding genealogical record collections and storing them in a spatial index (browse our blog if you want to know more). Searches are based on a shape. We return collections which intersect or are contained within the search area.

During the early states of development, we used geonames to geocode collections to a point. This works well for collections which are associated with lower administrative levels such as cities, towns, and villages. However it really breaks down when you get the county, state, and country level.

The 1940 US Census is associated with the United States and would be assigned a point in Northern Kansas. Any queries within the US that aren't near that point won't return the 1940 US Census.

To solve this we need to geocode collections with a shape instead of a point.

OSM

OpenStreetMaps has the data we need, but it's extremely difficult to extract. The administrative hierarchy is not explicitly stored. Nominatim is used to solve this problem for OSM a Nominatim search only returns features. So a query for Knighton on Teme returns two bus stops but not the administrative boundary relation.

The Overpass API looked promising but it can't do fuzzy string matches. Overpass can only do exact or regex matches. We could use Overpass if there was an easy way to standardize place names. In other words, if OSM provided a way for us to standardize "Knighton on Teme, Worcestershire, England" to "Knighton on Teme CP, Malvern Hills, Worcestershire, West Midlands, England, United Kingdom" according to the OSM hierarchy then fuzzy string matching wouldn't be necessary.

Summary

What we need is a service which allows us to perform fuzzy string searches for a place (or administrative level) and retrieve it's boundaries.

We recognize that it will be difficult to obtain boundary data for the entire world. Thankfully we probably won't need to anytime soon. We only need data for areas of the world where genealogical records exist and genealogists do research.

It's looking like we will need to build our own service which indexes OSM in such a way that enables us to query for administrative boundaries. But we would really prefer not to. Is there any other way we can retrieve this data with existing services?


texts: What simple and light Java library to use for WGS84 to UTM conversion?

texts:
Geotools and Geotoolkit are huge and complicated.

CTS (1.3.3 from Maven) failed on me by not even being able to run its own test example (slightly modified https://github.com/orbisgis/cts/blob/master/src/test/java/org/cts/op/EPSGTransform.java gave me "This registry EPSG is not supported" while I did not even touch those strings).

So, are there any alternatives that are easy to use, not require additional libraries if possible (even loggers are annoying enough) and are free software?

All I want to do is turn some WGS84 coordinates to UTM32N (EPSG:25832) in Java.


texts: Classify LAS By Height tool is not yielding any results?

texts:
The classify LAS by height tool (CLHT) in ArcGIS says that it will only reclassify LAS points with class codes of 0 or 1 (http://pro.arcgis.com/en/pro-app/tool-reference/3d-analyst/classify-las-by-height.htm). I have used the Change LAS Class Codes tool to reclassify class 3, 4, and 5 into 0. I then use the CLHT to enter in new class codes (>31 class codes) and respective heights. When the tool finishes, there are no new classes added to my LASD.

Where am I going wrong?

arcpy.ClassifyLasByHeight_3d(lasd, "GROUND", "32 1;33 2;34 4", "NONE", "false", "DEFAULT", "false")


texts:
I have a set of data in the column geometry(geom) as WKB.

I am using ST_AsGeoJSON(ST_Transform(geom, 32633)) to transform WKB into meters.

my results are like this:

{"type":"LineString","coordinates":[[100,100,100],[300,300,300]]

these basically are x, y and z coordinates accordingly. Namely i am getting results with 3 numbers for each set

What i want to do is deduct from each set of numbers the first set. That would mean that i want to do -[100,100,100] from all the next sets, because the first one should be my position on a GPS and it should be 0,0,0.


texts:
Imagine the following array:

a = [10768,10764,10893,10763,10723,10792,0,0,0,0]

Using a function like "MEAN" in IDL, is there a way to calculate the mean value, ignoring the zeros? I know this can easily be done in a for-loop, but I would prefer a ready made function for the sake of simplicity and performance.


texts: GIS Database Migration Project (Oracle to ArcGIS)

texts:
I have been tasked with converting an Oracle Database into an ESRI file geodatabase.  The data represents a gas network.  The client's Oracle database has over 400 tables, 120 of which are spatial.  I'm going to use FME to convert the data into an ESRI File Geodatabase.  The client has supplied me with the geodatabase which contains feature datasets and feature classes which are empty (schema only).  The geodatabase is a modified version of the ArcGIS gas data model (pdf).  Essentially it contains the main feature dataset with a geometric network as well as other data the client is already using in ArcFM.

The problem I have is that the Oracle database contains many more tables than the ESRI schema.  There are tables in the Oracle database which have been created which do not have an equivalent table in the geodatabase.  Furthermore, for those tables that I can match, there are fields that do not match.  For example, in the original database, there may be a table which contains fields such as "street number", "street name", "city", etc and the destination database does not contain these.  Furthermore, in the valve table in the original database, there is a field called "Direction To Close" which contains a "C" for Counterclockwise and an "A" for Anticlockwise and some blank values.  In the destination database there is an equivalent field called "ClockwiseToClose" which has a domain set for "Yes/No" values.

That is just one field in one table.  And as I said, there are over 400 tables.  I'll need to determine all of the issues which may arise for all fields in all tables.  Obviously, there is a lot of work involved in setting this all up.

Does anyone have any experience in a similar project?  Can anyone give a realistic estimate of how long a project like this might take?  I'm thinking several months.

Obviously the domains and subtypes in the destination geodatabase will likely cause issues and the data will have to be cleansed in the process.  Also, there may be fields which have numeric values in one database and text/categorical in another (with the domains).  That is, what might be "A", "B", "C" in the original might have a domain which restricts the values to "1", "2", "3", for example.

What other traps might I come across in a data migration project like this?


texts:
I have met a small yet puzzling problem about saving and reading Spatial Pixel Data Frame.
I created a raster with over 10 million grids (resolution = 25000).
adm0 <- readOGR(dsn = path.expand("C:/Users/Abel/Desktop/district"),
                layer = "administrative0")

proj.crs <- proj4string(CRS("+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

adm0.newproj <- spTransform(adm0, proj.crs)

grid_2.5km <- raster(adm0.newproj, res = 2500)  # 13143792 grids

In order to add data into the raster, I further transformed this raster into SpatialPixelsDataFrame (I did not use SpatialPolygons because my computer is incapable of dealing with vector of this size).
grid_2.5km[]<- 1:ncell(grid_2.5km)

grid_2.5km.sp <- as(grid_2.5km, "SpatialPolygonsDataFrame")

At first I used writeOGR to save this grid file and readOGR to read it. But turns out by doing so all grids would be transformed to SpatialPoint, which is definitely not what I want.
I also tried writeGDAL. But I have no idea about the set the filename, drivername, type and many other parameters.
writeGDAL(grid_2.5km.sp, "C:/Users/Abel/Desktop/divide/grid_2.5km.tif",
          drivername = "Gtiff", type = "Float32" )

How can I save this SpatialPixelsDataFrame correctly and keep its original object, and how should I read this file, letting the opened file still having the object as SpatialPixelsDataFrame?


texts:
I'm currently working on S1-images, and I want to process PSI on them. I use ENVI 5.5 with SarScape.

Here are the steps I followed :

I opened 23+1 images (Import Data) and I cropped them using SAR sampling function. I used orbit files.
Step1 - Connection Graph, everything seems OK.
Step2 - Interferogram Processing. There I use a DEM downloaded on USGS website (SRTM 35sec).

Unfortunately, I got the message error : "Memory BAD_ALLOC (-1.0)". On my computer, I have 1To space, and 16Go RAM. I'm on Windows 7, 64bits.
I tried many things to solve this problem :

Modify the AOI, down to pwr files measuring less than 8MB
Modify the DEM (using the DEM provided by SNAP)
Select less images (just to check) and the same error appears.

I did not find any solution to solve this problem. I know this seems like a memory problem, but even with very small images I got the issue. I did not find any piece of solution.

Did someone already have this problem ? Have you any ideas on how I can solve this ?


texts: How to get WKT from a QgsGeometry object?

texts:
How can I get the wkt from a geometry object in python?
I'd like to be able to do something like the following:

from qgis.core import *
import qgis.utils
layer = qgis.utils.iface.activeLayer()
for feat in layer.getFeatures():
    geom = feat.geometry()
    #print something meaningful about the geom object

I've checked dir(geom), but there doesn't seem to be an option to get a wkt.


texts: Calculate fields referring to other layer?

texts:
I have to analyze different pathway segments based on a few variables:

Number of properties crossed (polygons)
Length of paths on public land (polygons)
Number of intersection with stream (lines)
Number of intersection with river (polygons)

Is there a way to calculate these information directly in the attribute table using the field calculator?


texts: How to reimport/copy a shapefile where attribute table can be subsequently edited without changing original shapefile in QGIS

texts:
I am having trouble when duplicating a shapefile and then going on to edit the attribute table of the duplicate manually. The changes I make are reflected in the original file - hence creating a duplicate doesn't create a separate shapefile. I have tried copying the shapefile manually in the folder but then I found the attributes weren't in the shapefile.

QGIS 2.8.9


texts: Joining attributes from point layer to line layer which the endpoint's geometry matches the point geometry

texts:
In QGIS 2.18 I have 2 layers, one is a multiline layer and the other is a point layer.
The point layer contain every endpoints of that multiline layer and saved some of the point's own attributes.
I'd like to create 2 new fields in the multiline layer which are "start point type" and "end point type"  and the type of start and end point of each line are recorded in the attribute table of point layer.
The question is that how to extract the content from other layer's attribute table(and with different geometry type).
I've installed the refFunction plugin and find some function may be useful like geomequals.
Following is something I've try in the field calculator but not working(take start point for example).
case
    when start_point($geometry) geomequals('LayerL','$geometry')
    then get_feature('LayerP', 'FieldP', '')
end

LayerL is the multiline layer, LayerP is the point layer and FieldP is the field of point layer I'd like to fill in the new field.


texts:
I am using dotspatial version 1.9 to reproject a shapefile form EPSG:2326 to EPSG:3857 and I am getting wrong results.
For a single point
x:835869.4538390070
y:815974.4885250400

the result is
x:12709407.9774734
y:2545678.93085044
when using proj4 string (and completely different when I am using the ESRI String)

The correct outcome (according to epsg.io online transformation tool) should be:
x: 12709681.01
y: 2545494.73

which is in accordance with what I get using PostGIS.

Any ideas why the dotspatial doesn't give me the correct result?

The code I am using is:

        Double x;
        Double y;
        x = 835869.4538390070;
        y = 815974.4885250400;

        String proj4_2326_str = "+proj=tmerc +lat_0=22.31213333333334 +lon_0=114.1785555555556 +k=1 +x_0=836694.05 +y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m +no_defs";
        String proj4_3857_str = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs";
        ProjectionInfo proj2326 = ProjectionInfo.FromProj4String(proj4_2326_str);
        ProjectionInfo proj3857 = ProjectionInfo.FromProj4String(proj4_3857_str);

        double[] xy = new double[2] { x, y };
        double[] z = new double[1] { 0 };

        Reproject.ReprojectPoints(xy, z, proj2326, proj3857, 0, 1);

        textBox7.Text = xy[0].ToString();
        textBox8.Text = xy[1].ToString();