texts: Importing multiple ASCII text files in ArcGIS as Rasters

texts:
I have a big number of ASCII text files which have the Esri ASCII Grid format but have been stored as .txt files rather than asc. I am trying to import them simultaneously as rasters to obtain mean values from the various files.

I have already tried to create a mosaic dataset and import the files as Rasters but I get a 999999 error (potentially because my files have a different extension? that is txt instead of asc.)

Please see below a sample of my dataset - that is one of the many ASCII text files that I need to import as raster


texts: Only my geometry layers are exported from spatialite to postgis

texts:
I have a spatialite database (v4.1.1 / SQlite v3.7.17) consisting of one geometry table and several attribute tables linked to the geometry table by foreign keys.  I want to export all tables to a postgis (v2.1.3) / postgresql database (v9.3).  Using ogr2ogr and the following code:

ogr2ogr -append -lco GEOMETRY_NAME=geom -lco SCHEMA=public -f PostgreSQL PG:"dbname=... host=... port=... user=... password=..." -a_srs EPSG:4269 c:\MyFile.sqlite

I am able to export the geometry table, but not the attribute tables.  Any idea why?  Is there an option of the ogr2ogr function I am missing?


texts: Tilemill and GeoTiff not working together

texts:
Loading in a geotiff with projection of EPSG:4326 isn't fully rendering:

The original:

The original in Tilemill:

When using gdalwarp to change the projection from EPSG:4326 to EPSG:3857, tilemill will
then load the full geotiff, but the geotiff has a black strip down the center.

The warped in Tilemill:

The command used for gdalwarp is:

gdalwarp -overwrite -s_srs EPSG:4326 -t_srs EPSG:3857 -r lanczos -srcnodata 0 -dstnodata 0 -of GTiff original.tif warped.tif

The question is, either how do we get gdalwarp to work correctly, without the black strip, or how do we get the original to show up correctly in Tilemill.


texts: Using NASA Above landcover imagery in place with Qgis

texts:
Landsat-derived Annual Dominant Land Cover Across ABoVE Core Domain is composite tiled classified imagery dataset covering the years 1984-2014. The format is elegant and clever but hard to use. Each tile file is a 31-band GeoTIFF, with each band providing one year of data for the period (band_1 = 1984, ...band_31 = 2014) and classified into 15 values (1-15 in integer).

Neither ArcGIS or Qgis can use the data very well out of the box as they only allow colour maps and/or attributes for single channel 8 bit rasters. We could explode the bands into distinct files, but for my area of interest that will yield 175*31=5,425 files. A management headache to say the least, and probably file-system performance limit as well not to mention the increase in storage needs.

My intuition says both Arc and Qgis (via gdal) have the necessary features to be able to effectively use these imagery in place via raster mosaic-dataset with function chains (arc) and virtual raster tables vrt with pixel functions (gdal, and arc to lesser degree). I have so far though not managed to figure out the necessary magic.
The end goal is to enable people to load a given year from the data source and have it displayed in the standard colour map. Ideally it also yield the corresponding code when interrogated in ArcMap/Pro or Qgis. It should also be straightforward to use an arbitraty year as input to a geoprocessing tool. All with the least amount of fuss and hairpulling. ;-)
I would prefer a single solution that can serve both ArcGIS and Qgis but as requested I'm splitting this into a separate question for each platform. This is the Qgis posting. The ArcGIS posting is here.
I'm going to self-answer with the best I've come up with so far for mosaic solutions, but they're just partial steps along the way and I'm not happy with them. Please feel free to improve or better.


texts: Using QGIS to make changes to OpenStreetMap?

texts:
I want to be able to use a familiar editor (QGIS) to save building footprints to OpenStreetMap as I edit them.

I have been building a layer in QGIS and marking the polygons that need to be uploaded. I don't want to use JOSM which feels like it was written in 1985, and I don't want to re-digitize each polygon in the ID editor, because I am 62 and would die of old age before that process was completed.

There is an ArcMap tool, but it seems like there should be something in QGIS.

My data sources are license compliant (imagery and LiDAR).


texts: Constructing R-tree for continuous data?

texts:
I want to construct an R-tree to index land-use data.

How do we index continuous data using an R-tree?


texts: Transforming file geodatabase feature classes from custom local coordinate system?

texts:
Background:  Started working for company that is in middle of system-wide IT upgrade including GIS data.  Existing in-house data is in a FileGDB, is in an assumed coordinate system and represents 1000 sq\mi of service area in about 150 feature classes. I do have very good common data between existing dataset and 'real world' data (PLSS section data mostly).

System:  Working with ArcGIS Desktop 10.0 with ArcInfo license level.

In old school command line ArcInfo I would have written an AML to run a Transform on the data coverages using the section corners for tic references.  Not sure what the process is for File GDB.  I looked into the Transform command, and it only appears to work at the Workstation software install level and not at the Desktop level, even though I have the Info license.  Not sure if the Transform would work on a GDB instead of a coverage in any case.

Is there an equivalent Transform command for File GDB featureclasses?
What are any other options for converting?

My dataset is currently not in any coordinate system - it's floating\assumed\in-house coords.  The data is good within itself but not at all related to any real world cood system.  The old ArcInfo Transform command would use ID matched tics to move\warp\stretch data into a new spatial location.  I do have data segregated in Feature Datasets with the File GDB.


texts: Cannot generate buffers in meters (even after "save as" in NAD83 / UTM 18N) using QGIS?

texts:
I am quite new to GIS and QGIS.

I have 800 points with lat/lon (some locations in Massachusetts) in a .csv file. I use AddLayer to import the CSV file and save it as NAD83 / UTM 18N projection (EPSG:26918). Then I close QGIS and reopen the new 26918 file with QGIS. I want to create buffer for the points, and I want 500 meters as the buffer radius. I use Vector -> Geoprocessing Tools -> Buffer, and I entered 500 in the field for "Buffer Distance" and saved the new shapefile. Then the following pictures shows... I wonder how can I make such a buffer with 500m radius for the points?

The problem is that the buffer is way larger than 500 m (probably 500 degrees). I have figure out how to deal with it: after adding the csv file as a layer, I should first saving it in the WGS84 projection. And then save the WGS84 layer as a NAD83 / UTM 18N layer. Then I can create 500 m radius buffer.


texts: Double Entries/ Errors with Digitizing Polygons

texts: I have a layer that has errors, i mean polygons have double entries (digitized twice), i can edit manually but there's a lot of attributes and it will take ages to clean the data. Is there any short cut key or tool that i can use to pick up these polygons.


texts:
I built a model using model builder which I want to run using a python script.

The error i get is:

  AttributeError: Object: Tool or environment  not
  found

Can anyone please help me solve this as I am new to python.

The script I found online is:

import sys, string, os, arcgisscripting
import arcpy
arcpy.ImportToolbox
arcpy.env.workspace = 'Path of my geodatabase'
arcpy.env.overwriteOutput = True
# Create the Geoprocessor object...
gp = arcgisscripting.create()

# Set the License...
gp.SetProduct("ArcView")

# Load required toolboxes...
gp.AddToolbox("Path of my toolbox")

# Set the Overwrite property to True
gp.OverwriteOutput = 1

# Run the model
# All arguments were supplied to model prior to saving
gp.MultiLine_Geocoder()


texts: Estimate total stream length, frequency and order in ArcMap

texts:
I am new to GIS, and I want to calculate the stream frequency and Strahler stream order in a watershed. I have already processed the DEM and there are 20 sub basins. A sample map what I want to get is attached.

Can someone give some guidance in this regard?

May be a link, steps or tutorial will be good.


texts:
I am trying to convert a process usually used in ArcGIS Tools to arcpy in a Jupyter notebook.  I have a very long Con statement that basically takes all the values of 10 rasters and adds them together.  I have tested a small snippet of the code to make sure the syntax is correct, but I am running into errors saying 'TypeError: IsNull() takes 1 positional argument but 3 were given'
# Import system modules
import arcpy
import os

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

from arcpy.sa import *

# Set overwrite option
arcpy.env.overwriteOutput = True

# Set the workspace for ListFeatureClasses
arcpy.env.workspace = r"G:\Projects\bdv_2020Q4.gdb"
print(arcpy.env.workspace)

# Set the extent environment using the Extent class.
arcpy.env.extent = arcpy.Extent(-180.0, 90.0, 180.0, -90.0)

# Set the cell size environment using a number.
arcpy.env.cellSize = 0.0083333333

#Con statement
outputRaster = Con(IsNull(arcpy.Raster("IUCN_1a"),0,"IUCN_1a")) + Con(IsNull(arcpy.Raster("IUCN_not_reported"),0,"IUCN_not_reported"))
outputRaster.save("comb_test2")


texts: ArcGIS combine labels/annotations by feature

texts:
I am trying to find a way to combine multiple data points into one based on the same latitude longitude in a dataset. Basically too many labels appear for one point.

Example data

Name     Latitude     Longitude     Time
John     12.34567     -12.34567     10:42 AM
Marley   98.76543     -98.76543     11:52 AM
Alex     98.76543     -98.76543     12:00 PM
Heather  12.34567     -12.34567     9:16 AM
Angie    12.34567     -12.34567     5:47 PM

I am trying to combine this is a way so that the time labels will appear as one for the same latitude and longitude points. Is there an easy automated way to accomplish this?

Example output

Name                        Latitude     Longitude     Time
John \n Heather \n Angie     12.34567     -12.34567    10:42 AM \n 9:16 AM \n 5:47 PM
Marley \n Alex               98.76543     -98.76543    11:52 AM \n 12:00 PM

It does not need to be in the above format with newlines. Just trying to think of ways to do this.


texts: map GeoPandas dataframe with folium no results

texts:
I am trying to display on a folium map a set of polygons from a geodataframe. The code returns no errors but the polygons do not appear.
I have checked the projections ('epsg:3857') and the coordinates in the geo-json seem to be ok.
the code used is as follows:
import folium
map = folium.Map(location = [22, 55], tiles='cartodbpositron' , zoom_start = 5)

draft = hot_spots_draft
draft.crs = {'init' : 'epsg:3857'}

for _, r in draft.iterrows():
    #without simplifying the representation of each borough, the map might not be displayed
    #sim_geo = gpd.GeoSeries(r['geometry'])
    sim_geo = gpd.GeoSeries(r['geometry']).simplify(tolerance=0.001)
    geo_j = sim_geo.to_json()
    #print(geo_j)
    geo_j = folium.GeoJson(data=geo_j,
                           style_function=lambda x: {'fillColor': 'orange'})
    folium.Popup(r['ID']).add_to(geo_j)
    geo_j.add_to(map)
map

To test the geodataframe I used the code below and I got a simple plot
hot_spots_draft.plot(figsize=(10, 50), alpha=0.5, edgecolor='k')

I also exported the geodataframe to a shapefile with no issues.


texts: Changing the Processing plugin icon - QGIS

texts:
I have built a processing plugin which shows up in the QGIS Processing Toolbox. However, I am struggling trying to change the icon. Usually I do:
pyrcc5 -o resources.py resources.qrc
and it updates perfectly.
However, I do not have a
resources.py
file so I cannot adopt this approach. Is there another way to update it?


texts:
Has anyone ever heard of any tools in either Arcmap, QGIS, Mapwindow GIS etc, that could create a line plan (Survey Lines) across a polygon (Survey Area)?

A form of input such as distance between lines and minimum length would hopefully be required. If this tool could be tailored to create the shortest, so most efficient line plan, then that would be even better. Obviously a form of automation rather than just manually offsetting poly-lines is required- as that can get very monotonous over large areas!

The image attached shows the kind of output I need, however the Survey Areas are usually much more irregular!


texts: How to convert KML with tracks to Shapefiles?

texts:
I downloaded a KML track from Google Latitude. The service allows to download the tracked route from the user's profile providing a link such as:

https://maps.google.com/locationhistory/b/0/kml?startTime=1318716000000&endTime=1318802400000

ogr2ogr

I would like to convert the .kml file into a Shapefile. I am running Ubuntu. I am a bit familiar with QGIS and ogr2ogr. Therefore, I tried the following command as usual:

$ ogr2ogr -f "ESRI Shapefile" example.shp example.kml
ERROR 4: No layers in KML file: example.kml.
FAILURE:
Unable to open datasource `example.kml' with the following drivers.
[...]

Quantum GIS

QGIS states "invalid data source" when I open the .kml file as a vector layer.

Kml2Shp Online

I also tried Kml2Shp Online without success. It is not able to read any information either:

Entities found:
# Points: 0
# Paths: 0
# Inner Polygons: 0
# Outer Polygons: 0

Google Earth / GPSBabel

Furthermore, I found out that Google Earth is not capable of exporting GPX. Google refers to GPSBabel which converts KML to GPX using the following command:

$ gpsbabel -i kml -f ~/Desktop/example.kml -o gpx -F ~/Desktop/example.gpx

However, in my case it outputs a GPX file without location data.

GPS Visualizer

As a first success I found the website GPS Visualizer. Uploading and converting the file their actually produces a valid GPX file.

Still my question remains to be open:
Is there a command line tool that is capable of converting the file? I do not want to upload the tracking data to yet another website for the conversion.


texts: Comparing GRASS v.Overlay/QGIS Intersect Speed intersecting line and polygon?

texts:
Im currently trying to convert my scripts from ArcGIS to QGIS but one aspect that is extremely slow is trying to intersect a polygon with a line feature. I feel as though there must be a better way to accomplish trying to clip the voronoi polylines with the original polygon itself i.e. :

by

These are some of the bench tests that i've done with these two layers.

ArcGIS 10.1 Intersect
2min 24sec

Grass 7 v.Overlay
quit after ~1hour (~70% Selecting lines...) - Currently not supported in processing (Sextante)

Grass 7 v.split + v.Overlay
~2min + 45min

QGIS 1.9
quit after ~3hours

ogr2ogr -clipsrc
quit after > 1hour

SAGA GIS Intersect Line-Polygon
~15min -- acceptable but still 7x slower than ArcGIS and produces wrong results? The images below show SAGA GIS Intersect(with dangles removed) and same area with ArcGIS Intersect - Clearly SAGA did not intersect the line and polygon correctly as there should be no connection of the lines. These errors are scattered throughout the new output!!! In addition Line-Polygon does not retain the polygon intersect attributes to the line unlike the other methods mentioned above...

I have placed the original dataset in a GitHub repository here so hopefully others can make their own benchmarks and thereby provide suggestions?

I notice that ArcGIS Desktop splits (cracking) the featureclass so I do not know if that is another possibility to improve intersection speed within one of the opensource alternatives mentioned above?

***Note that the reason I mention QGIS here is that I have made some scripts within the processing toolbox (formally known as Sextante) and i'd like to remain in that environment if at all possible.

System - Intel Core i7-2760QM CPU @2.40GHz and 8GB ram


texts: Placing several atlas' features into a single page

texts:
All the documentations and tutorials I've been reading about QGIS' Atlas are focused on creating an atlas in which every page is a feature. Or in other words: where every feature is placed in its own single page.

I am wondering if it could be possible to create an atlas in which all their features were placed within the same page, in order to accomplish something like this:

This image has been taken from this magnificent R tutorial and displays population evolution in an array of maps within the same page. I thought that could be achieved in QGIS by creating an atlas where each year is a different feature, but Atlas' default behaviour would place each year in a different page, which is not exactly what I want, as I am sacrificing the power of seeing an evolution in just an image.

Any clues on that?


texts: How to create more than one mean coordinate in QGIS?

texts:
I have a point shapefile with several thousands of points in it.

I've split the points up into spatial groups using a vectorgrid and pointcount.

I wish to create a mean coordinate for each group of points. However, where there are more than 200 points in each group, I wish to add two new points that are evenly distributed amongst the group.

Is it possible to create two or more new points using the 'mean coordinates' tool in QGIS?

Further issues i'm having with this method:

How do I get my new point (generated by the mean coordinate) to 'snap' ontop of a pre-existing point
Groups are created based on their pointcount. Where I have two groups of points that are spatially apart, yet with the same pointcount, each point is placed in the same group - is there a way to alter this - perhaps with a weighted mean?

Any suggestions are greatly appreciated, thanks.


texts:
In Styling GeoJSON by Attributes with Cesium I asked about styling a GeoJSON in Cesium. Now that i have resolved that, I am trying to animate my Cesium scene. I have added a clock to my Cesium scene with the relevant range.

Here is the code :

var clock = new Cesium.Clock({
            startTime: Cesium.JulianDate.fromIso8601('1979-01-01'),
            stopTime: Cesium.JulianDate.fromIso8601('2014-01-01'),
            clockRange: Cesium.ClockRange.LOOP_STOP,
            clockStep: Cesium.ClockStep.SYSTEM_CLOCK_MULTIPLIER,
            multiplier: 7889238,
            shouldAnimate: true
        });
        var viewer = new Cesium.Viewer("cesiumContainer", {
            imageryProvider: new Cesium.IonImageryProvider({assetId: 3845}),
            clockViewModel: new Cesium.ClockViewModel(clock),
            shouldAnimate : true
        });
        var source1 = 'part of path';
        var source2 = 'part 2 of path';
        var year_time = 'year from current clock time';

        getLayer(year_time);

        function getLayer(year) {
            var promise = Cesium.GeoJsonDataSource.load(source1 + year + source2);

What I need is to get the current time value from the clock as a year and pass that to a function that creates and styles the GeoJSON.

Does anyone know if this is possible?

UPDATE:

I have been able to find a way to get the year as the time changes but for some reason when i try to pass it along to the getLayer function it get an error.
"SyntaxError: Unexpected token < in JSON at position 0"
Here is my updated code sample:

            var clock = new Cesium.Clock({
            startTime: Cesium.JulianDate.fromIso8601('1979-01-01'),

            currentTime: Cesium.JulianDate.fromIso8601('1979-01-01'),

            stopTime: Cesium.JulianDate.fromIso8601('2014-01-01'),
            clockRange: Cesium.ClockRange.LOOP_STOP,
            clockStep: Cesium.ClockStep.SYSTEM_CLOCK_MULTIPLIER,

            multiplier: 15778476,

            shouldAnimate: true
        });
        var viewer = new Cesium.Viewer("cesiumContainer", {
            imageryProvider: new Cesium.IonImageryProvider({assetId: 3845}),
            clockViewModel: new Cesium.ClockViewModel(clock),

            shouldAnimate: true,
            skyAtmosphere: false,
            skyBox: false

        });
        viewer.camera.flyTo({
            destination: Cesium.Cartesian3.fromDegrees(-56.5, 45.45, 450000.0),
            orientation: {
                pitch: Cesium.Math.toRadians(-50.0)
            }
        });

        var source1 = 'GeoJSON path part 1';
        var source2 = 'GeoJSON path part 2';

        var year_time = '';
        viewer.clock.onTick.addEventListener(function (clock) {
            year_time = Cesium.JulianDate.toIso8601(viewer.clock.currentTime).substr(0, 4);
        });
        console.log(year_time);
        getLayer(year_time);

        function getLayer(year) {
            var promise = Cesium.GeoJsonDataSource.load(source1 + year + source2);

Now I think I'm getting close but I believe that I'm going to need a way to only call the layer and maybe remove the past layer. I'm pretty stuck right now so any advice would be welcome!


texts: Getting/downloading shapefile from OpenLayers/Web page and saving to local disk

texts:
I have developed a web application where i am using GeoServer and OpenLayers. I have displayed some Layers over Google map in a ASP.Net Web Page. I want to know is there any Method or Script available in OpenLayers for downloading and saving shapefile (the shapefile which is displayed over Google map in Web page) from web interface.

Actually the user wants to download the shapefile from Web interface and open it in his Dektop GIS (ArcView) system.


texts: Improving performance of Eliminate tool in ArcGIS Pro?

texts:
My data is very large (polygon feature class with 884187 features), thats why Python geoprocessing tools (which use geoprocesses like Intersect, Buffer etc.) don't work due to performance issues.

So my idea was to use the "Eliminate" tool in ArcGIS Pro to reduce the amount of sliver polygons. Problem is now, that this tool is also not working because the data is too large (it always stucks while "Reading feature ..."). I know that e.g. Intersect uses tiled geoprocessing automatically. Is there a way to do this manually (like a tool that splits the data into smaller parts), or is there a better way to remove sliver polygons, than "Eliminate"?

My data is a landuse map with many single part polygons:

Technical Data:

PC: 64bit Windows 7, 8gb ram, 2,8GHz intel i5
Hard Drive(on which data is safed): HGST Travelstar Z5K500 HTS545032A7E380,   320 GB, 2.5", SATA 3Gb/s - 5400 rpm - Buffer: 8 MB, Main Partition
Software: ArcGIS Pro 2.1.0
Data Source: Type: Shapefile, Projection (ETRS 1989 LAEA), Resolution (minY: -9.526.700m maxY: 15.968.500m minX: -8.462.600m maxX: 17.068.600m)


texts: What is the lowest and highest possible pseudo-range?

texts:
I am trying to compact the required data that needs to be sent to the server to post-process the receiver's position. To do this, I figured it'd be possible to offset the pseudo-range by a fixed value. This would decrease the actual value that needs to be transmitted. The idea behind this is that the satellite is always a few thousand kilometers away, no matter where it is in the world.

The pseudo-range is transmitted as a double, which is 64 bits long. For my application, every byte I don't need to send should be scraped. It seems to me that a 5-byte unsigned value should be enough to store the pseudo-range offset. It allows me to adjust the pseudo-range by 1100 km in mm precision. Possibly a 32-bit floating point value would suffice, even.

I have the freedom to use my own format for storing the pseudo-range value, as long as it requires less than 8 bytes.

Therefore I need to know, what is the range of the pseudo-range values one can receive?


texts: PostGIS - Optimize Point Outside polygon query for 4 million points

texts:
I am trying to do a 'point outside polygon' query for GBR with ~ 4 million points.
Most of these points are inside the boundary but a handful are falling way outside the boundary.

The country shape boundary (A0 level) has been downloaded from
http://www.diva-gis.org/gdata/
Point data is a set of Latitude and Longitude

Below steps are done to pre-process point data -

Add primary key GID (ALTER TABLE point_gbr ADD COLUMN gid serial PRIMARY KEY)
Add GEOM (ALTER TABLE point_gbr ADD COLUMN the_geom geography(POINT,4326))
Update Geom (UPDATE point_gbr SET the_geom = ST_SetSRID(ST_MakePoint(longitude,latitude),4326))
Create Spatial index (CREATE INDEX pt_idx_geom ON point_gbr USING GIST(the_geom))

For Polygon data -

Uploaded the A0 shapefile using 'Postgis3 Shapefile and DBF loader' and set the SRID to 4326.
While upload I also create an index on GEOM column.

I executed the below query to find the points outside the boundary only if they fall beyond 500m, but it is still running for 6+ hours. Not sure what am I doing wrong?

SELECT
  ST_Distance(a.the_geom::geography, b.geom::geography) as distance_metre, a.*
FROM
  point_gbr  a , polygon_gbr  b
WHERE
  a.gid  NOT IN
    (SELECT a.gid
     FROM point_gbr  a ,  polygon_gbr  b
     WHERE ST_Within(a.the_geom::geometry, b.geom)
    )
 AND ST_Distance(a.the_geom::geography, b.geom::geography) > '500'

The laptop has a RAM of 24 GB and 4 cores and is running PostgreSQL 11 with PostGIS 2.5.

Also, is there a better way to do this for similar complex countries like Russia, Mexico?


texts: Checking accuracy of supervised classification method in ArcGIS Pro

texts:
I have a question about my image classification method in ArcGIS Pro 2.4.

I'm looking at a field of grain sorghum and doing a project that requires me to define certain features within the image for a conditional extraction. I began with this (multispectral) raster:

Using ArcGIS Pro's image classification wizard, I used training samples to classify the image into 5 different classes: Soils, shadows, Leaves, grain heads, and ground targets (coordinates are blanked out).

After computing the NDVI (and other indices), I then used a conditional statement to make a third raster to only include the "leaves" class, which is the only feature in the image I want to use.

How can I check the statistical accuracy of this method?

I've been looking into this, but being new to image classification, I'm not quite sure how to do this. This would be helpful so that I can compare it to other methods.


texts: Which county boundaries does Excel Power Map support?

texts:
I've got a bunch of UK data that is already aggregated to "county" level, that would merrily work in Mappoint 2006, however the borders it seems to use appear to be "odd", and some places just don't seem to be recognised.

What I don't understand is that we appear to have a boundary for the West Midlands (which it's matched against Walsall), but not Avon (which fell out of existence at the same time as the West Midlands). Lincolnshire appears to have merged the portions that were in Humberside, Rutland has popped into existence, and most of South Wales doesn't match up.

If there a way of influencing which sets of county borders to use, a way of making it only use borders that are in use at the same time, on the same map, or at least a list of the supported county names?


texts:
I can successfully create this featureclass in the ArcGIS 10 default geodatabase and add it to the project, but it can't be edited.  I have tried to add a feature using FeatureClass.CreateFeatureBuffer() but get the old cryptic: {"Error HRESULT E_FAIL has been returned from a call to a COM component."}

If I try to add a feature with the editing toolbar in ArcMap:

  The create feature task could not be completed.

The featureclass DOES have a spatial reference that matches the data frame and the other layers in the project.
Am I missing something here?  Should be easy, any ideas?

private IFeatureClass MakePointFC(IFeatureLayer pLineLayer)
{
    try
    {
        IGxCatalogDefaultDatabase Defaultgdb = ArcMap.Application as IGxCatalogDefaultDatabase;
        Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory");
        IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);
        IWorkspace pWorkspace = workspaceFactory.OpenFromFile(Defaultgdb.DefaultDatabaseName.PathName, 0);

        IFeatureWorkspace workspace = pWorkspace as IFeatureWorkspace;
        UID CLSID = new UID();
        CLSID.Value = "esriGeoDatabase.Object";

        IFields pFields = new FieldsClass();
        IFieldsEdit pFieldsEdit = pFields as IFieldsEdit;
        pFieldsEdit.FieldCount_2 = 2;

        IGeoDataset geoDataset = pLineLayer as IGeoDataset;

        IGeometryDef pGeomDef = new GeometryDef();
        IGeometryDefEdit pGeomDefEdit = pGeomDef as IGeometryDefEdit;
        pGeomDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPoint;
        pGeomDefEdit.SpatialReference_2 = geoDataset.SpatialReference;

        IField pField;
        IFieldEdit pFieldEdit;

        //pField = new FieldClass();
        //pFieldEdit = pField as IFieldEdit;
        //pFieldEdit.AliasName_2 = "ObjectID";
        //pFieldEdit.Name_2 = "ObjectID";
        //pFieldEdit.Type_2 = esriFieldType.esriFieldTypeOID;
        //pFieldsEdit.set_Field(0, pFieldEdit);

        pField = new FieldClass();
        pFieldEdit = pField as IFieldEdit;
        pFieldEdit.AliasName_2 = "SHAPE";
        pFieldEdit.Name_2 = "SHAPE";
        pFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;
        pFieldEdit.GeometryDef_2 = pGeomDef;
        pFieldsEdit.set_Field(0, pFieldEdit);

        pField = new FieldClass();
        pFieldEdit = pField as IFieldEdit;
        pFieldEdit.AliasName_2 = "Distance";
        pFieldEdit.Name_2 = "Distance";
        pFieldEdit.Type_2 = esriFieldType.esriFieldTypeDouble ;
        pFieldsEdit.set_Field(1, pFieldEdit);

        string strFCName = System.IO.Path.GetFileNameWithoutExtension(System.IO.Path.GetRandomFileName());

        IFeatureClass pFeatureClass = workspace.CreateFeatureClass(strFCName, pFieldsEdit, CLSID, null, esriFeatureType.esriFTSimple, "SHAPE", "");
        return pFeatureClass;

    }
    catch (Exception ex)
    {
        return null;
    }

}


texts: Collapse dual lines to centerline doesn't exist in ArcGIS Pro?

texts:
I have ArcGIS Pro build 3308 (Ver. 1.1) and wish to collapse two parallel lines into one, however it appears that this function no longer exists - I have also tried using python for this and received an error stating

AttributeError: 'module' object has no attribute 'CollapseDualLinesToCentreLine_arc'.

Was the function that allowed this to be done taken out in the pro version of ArcGIS?


texts: PostGIS latitude/longitude vs longitude/latitude confusion

texts:
I have been importing polygon data on PostGIS but started to wonder if I had input latitude/longitude in the right order. On PostGIS doc there is a page https://postgis.net/2013/08/18/tip_lon_lat/ that clearly states spatial database records longitude as x and latitude as y, which after reflection does make sense. It's also written on https://postgis.net/docs/ST_MakePoint.html "For geodetic coordinates, X is longitude and Y is latitude".
My confusion comes from the fact that I did import data with latitude/longitude as x/y. For example:
Text: POINT(36.060409 -79.816702)
X: 36.060409
Y: -79.816702

a point near Greensboro in the US.
So according to PostGIS, X is longitude and Y is latitude, but if I plot this it gives me a point in Antarctica. If I plot on Google map, as is (.i.e 36.060409, -79.816702), which does get coordinates as latitude,longitude, then the location is correct.
I am confused about how PostGIS records the data and how am I supposed to input it?
I use PostgreSQL 12, PostGIS 3.1.0 and SRID 4326


texts: Proper output while using geom.buffer in PyGIS?

texts:
I am developing a plugin in QGIS using PyQt5. I have used the following code to create a buffer of different distance for different features in a layer.

layers = [layer for layer in QgsProject.instance().mapLayers().values()]
drainIndex = self.dlg.comboBox_4.currentIndex()
drain = layers[drainIndex]
path=os.path.dirname(unicode(drain.dataProvider().dataSourceUri()))
crs=drain.crs()
outLayer = QgsVectorLayer(path,'Buffer1',"ogr")
out=os.path.join(path,'Buffer1.shp')
prov = outLayer.dataProvider()
outLayer.startEditing()
   for feat in drain.getFeatures():
        inAttr = feat.attributes() # Input attributes
        inGeom = feat.geometry() # Input geometry
        Segments=20
        bf_inGeom = inGeom.buffer(feat['B1'], Segments)
        poly=bf_inGeom.asPolygon()
        outGeom = QgsFeature()
        outGeom.setGeometry(QgsGeometry.fromPolygonXY(poly)) # Output geometry
        outGeom.setAttributes(inAttr) # Output attributes
        prov.addFeatures([outGeom])
        QgsVectorFileWriter.writeAsVectorFormat(outLayer,out,'utf-8',destCRS=crs,driverName='ESRI Shapefile')
    outLayer.commitChanges()

It is generating buffer correctly around the polygon, which is an input. But also it gives some extra features which are the buffers on the input polygon.

Therefore it is giving double features than expected. One is the buffer outside input polygon and another is buffer on the input polygon. It is shown on the image given below.

As shown on the picture green coloured polygon is an input file. Red lines are buffers. I want buffer which is outside of input file.

What changes are recommended in the code?


texts: HTML pop-up in ArcGIS 10.1 by using HTML code within a column

texts:
I have a layer in ArcGIS 10.1 having a column in the attribute table filled in with HTML code (<p class=MsoNormal> A number of similar..), I was wondering whether it is possible to display a pop-up message based on that HTML code.

For example, by using the HTML Pop-up tool it is the closest I've been, but how do I tell ArcGIS that a specific column should be read as HTML code and not as simple text?.


texts: Exporting multiple maps on rendercomplete in OpenLayers

texts:
I have this project where I have 4 different maps next to each other. These maps have different layers on top of them so they can be compared. Now I want the user to be able to export those 4 maps to a pdf file. I found some documentation on how to export a map to pdf: Documentation.
I edited this down to make it compatible in my react project. At this point I am able to export one map. But not all 4 maps in one pdf file. This is because a map can only be added to a pdf once the map is rendered. But I can only check the "rendercomplete" on one map at a time. Is there a way to check if all maps in a Array are rendered and then execute the code to export to a pdf? In my example i have the this.Maps[0] just checking the first map in the Array. I need to check all maps in that Array and when they are all rendered the code needs to be executed.
  exportMap() {
    this.Maps[0].getMap().once("rendercomplete", () => {
      var pdf = new jsPDF("landscape", undefined, "a4");
      var viewResolution = this.Maps[0].getMap().getView().getResolution();
      var mapCanvas = document.createElement("canvas");
      var size = this.Maps[0].getMap().getSize();
      mapCanvas.width = size[0];
      mapCanvas.height = size[1];
      var mapContext = mapCanvas.getContext("2d");

      this.Maps[0]
        .getMap()
        .getViewport()
        .querySelectorAll(".ol-layer canvas")
        .forEach((canvas) => {
          if (canvas !== undefined) {
            if (canvas.width > 0) {
              console.log(canvas);
              var opacity = canvas.parentNode.style.opacity;
              mapContext.globalAlpha = opacity === "" ? 1 : Number(opacity);
              var transform = canvas.style.transform;
              // Get the transform parameters from the style's transform matrix
              var matrix = transform
                .match(/^matrix\(([^]*)\)$/)[1]
                .split(",")
                .map(Number);
              // Apply the transform to the export map context
              CanvasRenderingContext2D.prototype.setTransform.apply(
                mapContext,
                matrix
              );
              mapContext.drawImage(canvas, 0, 0);
            }
          }
        });
      pdf.addImage(mapCanvas.toDataURL("image/jpeg"), "JPEG", 0, 0, 297, 210);
      // Reset original map size
      this.Maps[0].getMap().setSize(size);
      this.Maps[0].getMap().getView().setResolution(viewResolution);
      document.body.style.cursor = "auto";
      pdf.save("map.pdf");
    });
    this.Maps[0].getMap().renderSync();
  }


texts: Installing exifread for QGIS photo2shape plugin for novice

texts:
I'm trying to install exifread so I can use the photo2shape QGIS plugin.

I know nothing about python and I haven't been able to get it to work using the answers from previous questions.

Can anyone please help with a set of clear step by step instructions for a novice?


texts: Geographic vs. projected raster cell size area

texts:
I am wondering the best way of knowing how to define the area of each raster cell in global data. Specifically, I downloaded all ~400 tiles of MODIS land cover data, mosaiced them into one, and now am summarizing land cover areas at admin levels. The data come downloaded as custom spheroid sinusoidal (~463m resolution). Each subset appears to contain cells all of the same size. My question is after I mosaic the data into one raster, if I changed the projection to WGS84, would I then have to account for the cell sizes changing with latitude? Does it depend on what the projection the original data were in? I can't think of any way to test this, since when looking at the data in ArcMap the cells all measure the same, using the measure tool, which may or may not just be a product of the dataframes coordinate system. Would converting a subset of the raster to shapefile and then calculating the area tell me anything?

Overall my question is a basic one, does every raster in a geographic coordinate system by its nature have cell sizes that change with latitude? Or are there exceptions?


texts:
Can raster maps be stored in a SQL Server spatial database?

And if so, how?

... I'm aware of VARBINARY and FILESTREAM for images, but can only see GEOMETRIC and GEOGRAPHY for vector data


texts: Specifying tile background color with Split Raster of ArcGIS Desktop?

texts:
In ArcMap Desktop 10.2.2, I have a Erdas Imagine, please see the picture, its NoData is displayed as blue.
I can use "Split Raster" tool(Opened by clicking: menu Geoprocessing >> ArcToolBox >> Data Management Tools >> Raster >> Raster Processing >> Split Raster)  to split this map into tiles.
however, in the splited tiles, the backgrond is black, (what I expect is background is blue).

How can I specify a color that nodata pixel will be rendered ?

I guess that the color may be related to the "NoData value" option when I runing "Split Raster", however, I searched all web but can not find any clues how to set this option.


texts:
We would like to improve our bike routing to also incorporate, for example, the information on the surface of the used streets. I don't quiet get how I should use the Osm2Po plugin in mechanism. I've set up my project in Eclipse and have a working development environment with a custom WayTagResolver (not sure if that is really needed) and a custom PostProcessor that in theory should write the additional information in the SQL file.

Is there a tutorial anywhere that can help me?


texts: Changing null value to specific numbers in ArcGIS for Desktop?

texts: I have a Excel file showing annual precipitation and I add it on arcmap. When I tried to check it's attribute table the field calculator is passive. I want to change all the Null value to specific number for analysis. How can I solve it please?


texts: How to open an Esri Personal Geodatabase (*.mdb;) in QGIS 3.4 of 64 bit in Windows 10 Pro?

texts: I am a learner. While I start to load Esri personal database, QGIS 3.4 does not support it.


texts: Given a lat/long point, determine if it's land or water?

texts:
I have a PostGIS-enabled PostgreSQL 9.1 installation with the 10M physical dataset from Natural Earth imported. It's been converted to 900913, and in QGIS I can visualize the data. So the plumbing's all set.

I have a very simple task: I'd like a SQL query to determine if a given point is land or water. Here's a very simple use case: user mouse-downs on a slippy map, the server takes that point and replies "that's land" or "that's not land." I feel that this is surely accomplished with the ne_10m_land layer (with enough accuracy for my purposes).

I'm a newbie to GIS, so while I've tried to do my homework I feel like I must be missing something really obvious to accomplish this...


texts: How could I compare polygons with appropriate circle?

texts: My problem is that I have a polygon shapefile with 1034 records which are the Greek municipalities. My goal is to compare the shape of every record with the circle. In other words, I want to observe the deviation of the polygons from the perfect shape-circle.


texts: Get single popup onclick different WMS layers using OpenLayers 3

texts:
I am trying to develop a web mapping application. I have more than two WMS layers and trying to GetFeatureInfo onclick the layer. I tried the answer posted on this link. But whenever I click on one layer three popups are displayed in spite of clicking only on one layer. I am stuck on how to resolve the issue.

Below is the code am working on:

    map.on('singleclick', function (evt1) {
    var coordinate = evt1.coordinate;
    content.innerHTML = '';
    var viewResolution =(view.getResolution());
    var url = '';
    content.innerHTML = '';
    layers.forEach(function (layer, i, layers) {
      if (layer.getVisible() && layer.get('name')!='Basemap') {
            url = layer.getSource().getGetFeatureInfoUrl(evt1.coordinate, viewResolution, 'EPSG:3857', {
                'INFO_FORMAT': 'text/html','FEATURE_COUNT': '300', 'WIDTH':500, 'HEIGHT':500,
            });

            if (url) {
              var left_f = "left_frame";
              content.innerHTML += '
You clicked here:
';
                overlay.setPosition(coordinate);
            }else{
                        content.innerHTML = '';
            }
        }
    });

});

Query:

GetFeatureInfo with a popup on click single layer or multiple layers. Is there any other way to solve this query?

I am new with OpenLayers and have basics of JavaScript. Below screenshots of output are attached.


texts: Inserting substring in ArcGIS Field Calculator

texts:
I'm trying to use the split tool in ArcGIS to add an underscore. Trying to make "ABC11111" look like "ABC_11111". Don't have much experience using Python or VB Script but this is what I have created with little success:

Split( [SURVSTR_ID] , "_", 3)(ABC).


texts: Iterate over GEE assets in Python script and mixing local and server side objects

texts:
I have below two questions that are related together:

1- Can I iterate over items in my GEE assets in Python script? If yes, how?

2- My understanding is that all the images and features (and their collections) are server side objects and I can not define an image based on a local file on my computer and then for example clip a Landsat image in GEE repository with it, am I right?


texts:
I want to display the image with its date as the layer name, but the following code does not work.
var s2 = ee.ImageCollection('COPERNICUS/S2_SR').filterBounds(roi).filterDate(startdate,enddate).first();
var date=s2.date()
Map.addLayer(s2, {}, ' '+date,true);


texts: What's the binary format of ENVI (Classic) ROIs?

texts:
I'm a Pythonist in an office full of ENVI/IDL users. I've created some tools to help bridge the gap, but one big stumbling block is ENVI Classic .roi files, which my colleagues use a lot.

I can easily import ENVI 5 ROIs into Python, since the text-based format is easy to parse. However, I don't have a good sense for the binary format of ENVI Classic .roi files, which make up the bulk of data interchange in my group. Does anyone have any code that reads or writes these files, or any insight on how one might go about it?

It's somewhat tedious to fire up ENVI 5 every time I want to go back and forth - I'd rather solve the problem at the source if possible.


texts:
I'm new to MapInfo and am struggling to import .gml OS mastermap files.

The .gml is created by converting the .gz data that OS provide to .gml using the following batch file:

FOR /R "D:\path\to\files" %f IN (*.) DO REN "%f" *.gml

I have tried downloading the OS schema's that they provide on their website which doesn’t change the error message that I get which is as follows:

Unable to access "C:\filepath\ss10ne.gml file

Can anyone recommend a solution/alternative?


texts: Renaming shapefiles from Cyrillic symbols to Latin ones using ArcPy

texts:
My goal is to rename shapefiles in 2 folders via transliteration using arcpy. That means that I want to replace Cyrillic (Russian) symbols with Latin ones. Here is my code, that doesn't work. What's wrong with it?

import arcpy

InWork =  ["C:\\Users\Vadim\Downloads\GML\esri", "C:\\Users\Vadim\Downloads\GML\stp"]

for ws in InWork:
    arcpy.env.workspace = ws
    datasets = arcpy.ListFeatureClasses()
    for fc in datasets:
        newName = fc.replace(u"абвгдеёзийклмнопрстуфхъыьэАБВГДЕЁЗИЙКЛМНОПРСТУФХЪЫЬЭ",u"abvgdeezijklmnoprstufh'y'eABVGDEEZIJKLMNOPRSTUFH'Y'E")
        arcpy.Rename_management(fc, newName[:-4])


texts: Accessing other layers using Arcade in ArcGIS Portal

texts:
I want to create a pop-up that uses multiple layers like in this blog post on Using FeatureSetBy functions in Arcade to drill-down to other layers and tables. However, the web map I am working with has many, many layers. Most of these are grouped into Group Layers and are part of a Map Image Layer published from an MXD.

Is there a way to use FeatureSetBy to drill into another layer within a group layer?

The globals to move up like $map or $datastore don't work like this and I am not sure what alternatives there may be.


texts: In Mapbox Studio how do I make an image for a marker be 400 pixels without getting clipped off?

texts:
For Mapbox Studio, I've prepared this Squid drawing I'm going to put in a body of water, similar to the ships in the 'Pirates Return' map.

I have sizes of the drawing in 100,200, and 400 pixels. However, the 400 pixel drawing gets cut off.  The last time I had this problem it was in the buffer settings. However, my buffer-size is now at the max.  So for an image applied to a marker, how do I make it overlap adjacent tiles?

Picture attached for clarification.


texts:
I've got 2 polygons layers (no gaps): IDUs and crops. What I want is to identify the 3 largest crops by % area for each IDU. This seems common enough but I couldn't find a way to do it all using tools, so I decided to try and learn basic arcpy to finish the job. I'm in ArcGIS 10.1 SP 1 (Advanced) and have done lots of PHP but this is my first try at python.

I found the Tablulate Intersections tool which works great and gives a table with all the overlapping crops with area and %area fields, but I want to add attributes to the IDU layer for the 3 largest crop codes - eg. LU_1=123, LU_2=201, LU_3=132. Here's the script I wrote. It works, but it's taking 100sec to process 160 test IDU records, and eventually it needs to process 100,000. I've got a 2013 Retina MacBookPro which flies through everything else I've thrown in in Arc so far, so something is obviously poorly done in my code.

import arcpy
idusFC = "c:/Envision/StudyAreas/PEI-Dev/NewIDUs/PEI.gdb/IDUs_Test"
cropTB = "c:/Envision/StudyAreas/PEI-Dev/NewIDUs/PEI.gdb/IDU_CROP2011_table_Sort1"
iduFields = ["CROP2011_A_Code", "CROP2011_B_Code", "CROP2011_C_Code", "IDU_ID"]
cropsFields = ["CROP_2011_GRIDCODE"]
print "Starting processing"
start = time.time()
with arcpy.da.UpdateCursor(idusFC, iduFields) as idus:
    for idu in idus:
        itr = 0
        iduID = idu[3]
        whereClause = arcpy.AddFieldDelimiters(cropTB, "IDU_ID") + " = " + str(iduID)
        with arcpy.da.SearchCursor(cropTB, cropsFields, where_clause=whereClause) as  iducrops:
            for iducrop in iducrops:
                idu[itr] = iducrop[0]
                itr += 1
                if itr > 2:
                    break
        idus.updateRow(idu)
end = time.time()
processTime = end - start
print "Processing Completed - {:.1f} seconds".format(processTime)

Is there another way I should do this so that it runs reasonably fast?


texts: Create point geometry in PostGIS from a X Y Z M coordinate format

texts:
I would like to get a PostGIS point table from a PostGreSQL table containing OpenDataKit data where the coordinates are contained into in a column in a X Y Z M format (M being I think the value of the GPS precision in meters). E.g:
10.384643333333335 -23.214940000000002 -0.8 5.0

My first question is: what should be the column type for my coordinates_column with such coordinates format? Just text?
Otherwise I try creating the geometry this way:
INSERT INTO schema.posgistable SELECT _id,column1,column2,...,ST_MakePoint(coordinates_column) as geom FROM schema.postgresql_odk_table

as ST_MakePoint seems to be the one able to retrieve X Y Z M data (but separated by ,)
but I got this error message:
ERROR:  parse error - invalid geometry

Any clue?


texts:
I'm using geoserver with geowebcache and I need to use the REST API to truncate tiles in a bounding box so they can be rebuilt. All the examples in the documentation are in curl and I'm using JQuery. I'd like to do the call with $.ajax. Here's what I'm trying:

var formData = new FormData();

            formData.append("threadCount", "01");
            formData.append("type", "truncate");
            formData.append("gridSetID", "EPSG:900913");
            formData.append("format", "image/png");
            formData.append("zoomStart", '00');
            formData.append("zoomStop", '15');
            formData.append("minX", "-12484306.954023438");
            formData.append("minY", '3991847.3646093756');
            formData.append("maxX", "-12445171.195546875");
            formData.append("maxY", '4030983.1230859384');

            var data = {
                "threadCount": "01" ,
                "type": "truncate",
                "gridSetID": "EPSG:900913",
                "format": "image/png",
                "zoomStart": '00',
                "zoomStop": '15',
                "minX": "-12484306.954023438",
                "minY": '3991847.3646093756',
                "maxX": "-12445171.195546875",
                "maxY": '4030983.1230859384'
            };

            var ajaxResponse = $.ajax({
                type: "get",
                url: "http://localhost/geoserver/gwc/rest/seed/MyWorkspace:cityJoin",
                contentType: false,
                data: data,
                processData: true,
                cache: false
            });

I have tried several variations of this code. I've tried sending the formData and the data, get and post, and I've changed the contentType variable between false, multipart-form, and json. I've also had processData set to true and false and not included at all. I've also tried appending all the data to the url and sending it that way. The data I'm using comes from a call the Geowebcache seeding page makes. I get one of three results based on the settings I use: 1. the call goes through but no tiles get truncated. 2. Illegal Invocation error. 3. Internal server error. I can't get it to run the truncate. Does anyone have an example of how to run the truncate with $.ajax or have an idea why my code isn't working?


texts:
I have a Python script that is working fine from Processing->Toolbox->User Scripts. However, I have put some user defined functions in "my_functions.py" that I can call from "Timing_Advance.py".

The problem is, when I edit the functions in "my_functions.py" in the QGIS script editor, I have to restart QGIS for the changes to take effect. This is perhaps because, when I restart QGIS, it makes a "my_fucntions.pyc" file and keeps it in the same location. The .pyc file is not update when I edit the user defined functions in "my_functions.py".

How can I update the my_functions.pyc file after editing, without having to restart QGIS?

UPDATE 24th Jan 2017
Looks like including the following lines isn't making any difference either

import sys
sys.path.append(r'C:\Users\EDEYARI\.qgis2\processing\scripts')
import py_compile
py_compile.compile(r'C:\Users\EDEYARI\.qgis2\processing\scripts\my_functions.py')

import my_functions

If I add some dummy lines like a "print" and re-run the script, the print doesn't give the desired output, till I re-start QGIS


texts:
Generate Near Table or Near can be used on a single Feature Class to find the nearest feature to each feature.

 arcpy.GenerateNearTable_analysis(featureClass, featureClass, pairsTable)

This results in two rows for each pair, when feature A is nearest feature B and when feature B is nearest feature A. Feature C may also be nearest feature B. I'm trying to come up with a table of unique pairs using arcpy (e.g. keep A-B, C-B, but remove B-A).

I started out doing this using a selection based on a cursor, and it seems to work but it is extremely slow (~20min for 1300 points), there must be a simpler way.

 arcpy.MakeFeatureLayer_management(pairsTable, "lyr")
 pairs_lst=[]
 with arcpy.da.SearchCursor(pairsTable, ["IN_FID", "NEAR_FID"]) as cursor:
      for row in cursor:
           pairs_lst.append([row[0],row[1]])
           if [row[1], row[0] not in pairs_lst:
                exp = '"IN_FID" = ' + str(row[0])
                arcpy.SelectLayerByAttribute_management("lyr", "ADD_TO_SELECTION", exp)
 arcpy.CopyFeatures("lyr", pairsTable2)


texts: I have an SBET (Smoothed Best Estimate of Trajectory) file that I would like to read in python.  I found these links that discuss options for reading/writing SBET files, but I was wondering if our community could offer some advice for handling SBET files in python?


texts: How to match name if layer name contain "\f"

texts:
One of my layer name is

"FullProbeWMS\fullprobewms[Heatmap] Trace Lines"

Where I check if this layer is present or not by using below code.

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.longName == "FullProbeWMS\fullprobewms\[Heatmap] Trace Lines":
        ispresnt = True
        islyr = lyr
        break
    else:
        ispresnt = False
        islyr = ""

However, when I run this code I always get values of "ispresnt" or "islyr" as "False" and "" because the name contains a character "\f".

How could I fix that?

I have tried many workarounds as changing values in string or single quotes etc.


texts: arcpy.mapping.ListLayers() without calling current mxd document

texts:
How can I use the arcpy.mapping.ListLayers() function without calling the current mxd? I'm hoping to run the Python script without having an ArcMap document open.

My setsymbology code runs extremely slowly as it is now. I could do it faster manually, and my computer is definitely not the limiting factor. I'm guessing it's slow because ArcMap is open. I'll show the code as it is now and the code that I'm trying to figure out a way for it to call an input directory/workspace rather than an mxd, but unfortunately it gets errors ('module' object has no attribute 'Listlayers').

# Import system modules
import arcpy
from arcpy import env

mxd = arcpy.mapping.MapDocument("CURRENT")

symbologyLayer = r"C:\VMshared\small_example_valley3\SnowColor2\snowdepthCOLOR.tif.lyr"

for lyr in arcpy.mapping.ListLayers(mxd):
    arcpy.ApplySymbologyFromLayer_management(lyr,symbologyLayer)


texts:
I'm trying to move a feature class on our network SDE that's made up of polylines as a CSV file to the local machine (or fileGDB but slow speed seems to be the theme regardless of output format)

import arcpy, csv
import datetime
from arcpy import env

def tableToCSV(input_tbl, csv_filepath):
    fld_list = arcpy.ListFields(input_tbl)
    fld_names = [fld.name for fld in fld_list]
    with open(csv_filepath, 'wb') as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(fld_names)
        with arcpy.da.SearchCursor(input_tbl, fld_names) as cursor:
            for row in cursor:
                writer.writerow(row)
        print csv_filepath + " CREATED"
    csv_file.close()

env.workspace = r"Q:\SDE\Direct Connection to delivery.sde"
desiredFC = arcpy.ListFeatureClasses("","","Topography")

for fc in desiredFC:
    if fc == 'Contours':
        print 'Starting to Process/' + fc
        startTime = datetime.datetime.now()
        out_csv = "C:\Users\CHOK\Desktop\Extract Data Method/" + fc[14:] + '.csv'
        tableToCSV(fc, out_csv)
        endTime = datetime.datetime.now()
        print 'Finished Processing/' + fc
        print 'It took ' + str(endTime - startTime) + ' to process ' + fc

raw_input('Press Enter to exit')

When I run the above script, I see that transferring 3.6 gb of data from the SDE to a CSV (or file geodatabase using arcpy copy methods) takes more than 3~4 hours to run.

I can see in task manager that the network speed utilised is usually between 0.3 Mbps to 2 Mbps but sometimes the speed jumps to 10~20 Mbps and other feature classes are transferred at around 5~7 Mbps.

Is there something I can try to improve the speed at which SDE stored on the network is copied over to a local drive or is this a lost cause due to the network currently in place?


texts: How to install xlw and xlrd libraries on MacOS so that QGIS Python can see them?

texts:
I experience a problem with XYtools and need to install Excel libraries for Python. I tried to install them directly from the terminal and it worked on the system, however, how to import them into QGIS?
This is how the error looks like:

import xlrd
Traceback (most recent call last):
File "<input>", line 1, in <module>
File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 572, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named xlrd

I tried to install them with pip:

python -m pip install xlwt==1.0.0
Requirement already satisfied (use --upgrade to upgrade): xlwt==1.0.0 in /Users/myname/anaconda/lib/python2.7/site-packages/xlwt-1.0.0-py2.7.egg

When I try to export file as xls from QGIS, the same error persists:

There is a problem with loading a python module which is needed to read/write Excel files. Please see documentation/help how to install python xlw and xlrd libraries.


texts: Joining attributes by location for multiple layers using PyQGIS

texts:
Using the script below, I've created a list of catchment area corresponding to a list of outlets. It works perfectly but I would like to keep the attributes of my outlets in my catchment area layers (output_3 in my script).
In order to do so I've added a new output (output4) to join each outlet layer to its corresponding catchment area.
import processing
import os # This is is needed in the pyqgis console also
from qgis.core import (QgsVectorLayer)
#provide path to outflow points where x, y are stored in fields named 'x', 'y' (or rename geom.asPoint() to match your coordinates fields)
path_outlet = 'F:/outlet.shp'

#provide path to filled and clipped DEM
path_filled_DEM = 'F:/mnt_surcreuse.sdat'

#provide path where to store temporary outputs from processing algorythims
path_to_temp_output1 = 'F:/upslope_raster{}.sdat'
path_to_temp_output2 = 'F:/polygonize{}.shp'
path_to_temp_output3 = 'F:/extract{}.shp'
path_to_temp_output4 = 'F:/final_output{}.shp'

#adds outflow point vector file to processing
outflow_points = QgsVectorLayer(path_outlet, "Outflow points layer", "ogr")

#iterates over all points in the outflow file. i is going to serve as polygone identifier
for i,f in enumerate(outflow_points.getFeatures()):

    #specifies geometry
    geom = f.geometry()

    #creates a raster file with upslope area
    output_1 = processing.run("saga:upslopearea", {'TARGET_PT_X': geom.asPoint().x(), 'TARGET_PT_Y': geom.asPoint().y(), 'ELEVATION': path_filled_DEM, 'METHOD': 0, 'AREA': path_to_temp_output1.format(str(i))})

    #converts raster to polygon
    output_2 = processing.run("gdal:polygonize", {'INPUT': output_1['AREA'], 'BAND': '1', 'OUTPUT': path_to_temp_output2.format(str(i))})

    #extracts only polygons with value 100
    output_3 = processing.run("qgis:extractbyattribute", {'INPUT': output_2['OUTPUT'], 'FIELD': 'DN', 'OPERATOR': '0', 'VALUE': '100', 'OUTPUT': path_to_temp_output3.format(str(i))})

    #add fields of each outlet layer to each corresponding upslope area
    output_4 = processing.run("qgis:joinattributesbylocation", {
        "TARGET": path_outlet,
        "JOIN": output_3,
        "PREDICATE": 'intersects',
        "SUMMARY": 0,
        "KEEP": 1,
        "OUTPUT": path_to_temp_output4.format(str(i))})

Unfortunately the script doesn't work and I get this error message instead:
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\script\ScriptEditorDialog.py", line 224, in runAlgorithm
    exec(self.editor.text(), _locals)
  File "<string>", line 50, in <module>
  File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\tools\general.py", line 106, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/PROGRA~1/QGIS3~1.10/apps/qgis/./python/plugins\processing\core\Processing.py", line 137, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.
Could not load source layer for INPUT: no value specified for parameter

(line 50 corresponds to the last line of the script)
Could someone help me solving this issue?
I am using QGIS 3.10.


texts: Using ArcMap Field Calculator to set Row 1 (Column D)= Row2 (Column B)?

texts: This is probably a fairly basic question but I have a long list of latitude and longitude values that represent telemetry points that will be the start point of the line.  Before I convert them to line, I need to add a column for the ending latitude and longitude values of the line.  This will be equal to the lat longs of the next point.  Each point is represented by a column "uniqueID"  This is one example of what I need done but there are several other instances in the database.  The example is shown in excel and is easily calculable with an "=B2" expression in excel but I need a formula for it in ModelBuilder.


texts: How can I optimize this ST_Distance query?

texts:
I am using ST_Distance in postgis to get the closest point in a table to each polygon in another table. This query is taking a long time to run (2 hours with 18,000 records in the polygon table and 7,500 records in the point table). I am calculating the distance in miles and also limiting the search to points that are within 100 miles (160934 meters) of each polygon. What else could I do to speed this up? The geometries are all in SRID 4326, and I have created spatial indexes on both tables. I'm running postgresql 9.5 with postgis 2.2.2. Here is the query:

SELECT DISTINCT ON(g1.gid) g1.*, g2.plant_name As closest_utility_name,
g2.primsource As closest_utility_type,
round(ST_Distance(g1.geom::geography,g2.geom::geography) * 0.000621371) As closest_utility_miles
INTO parcels_utility
FROM parcels As g1, powerplants_us_201603 As g2
WHERE ST_DWithin(g1.geom, g2.geom, 160934)
ORDER BY g1.gid, closest_utility_miles;

Explain provides this query plan:

"Unique  (cost=11006.38..11006.60 rows=45 width=713)"
"  ->  Sort  (cost=11006.38..11006.49 rows=45 width=713)"
"        Sort Key: g1.gid, (round((_st_distance((g1.geom)::geography, (g2.geom)::geography, '0'::double precision, true) * '0.000621371'::double precision)))"
"        ->  Nested Loop  (cost=0.15..11005.14 rows=45 width=713)"
"              ->  Seq Scan on parcels g1  (cost=0.00..1692.01 rows=18001 width=653)"
"              ->  Index Scan using powerplants_us_201603_gix on powerplants_us_201603 g2  (cost=0.15..0.51 rows=1 width=60)"
"                    Index Cond: (geom && st_expand(g1.geom, '160934'::double precision))"
"                    Filter: ((g1.geom && st_expand(geom, '160934'::double precision)) AND _st_dwithin(g1.geom, geom, '160934'::double precision))"

Based on the comment below I've updated the ST_DWithin function to use the geography type:

SELECT DISTINCT ON(g1.gid) g1.*, g2.plant_name As closest_utility_name,
g2.primsource As closest_utility_type,
round(ST_Distance(g1.geom::geography,g2.geom::geography) * 0.000621371) As closest_utility_miles
INTO parcels_utility
FROM parcels As g1, powerplants_us_201603 As g2
WHERE ST_DWithin(g1.geom::geography, g2.geom::geography, 160934)
ORDER BY g1.gid, closest_utility_miles;

It has been running now for about an hour, the query plan is below:

"Unique  (cost=612038648.76..612038650.68 rows=200 width=455)"
"  ->  Sort  (cost=612038648.76..612038649.72 rows=383 width=455)"
"        Sort Key: g1.gid, (round((_st_distance((g1.geom)::geography, (g2.geom)::geography, '0'::double precision, true) * '0.000621371'::double precision)))"
"        ->  Nested Loop  (cost=0.00..612038632.33 rows=383 width=455)"
"              Join Filter: (((g1.geom)::geography && _st_expand((g2.geom)::geography, '160934'::double precision)) AND ((g2.geom)::geography && _st_expand((g1.geom)::geography, '160934'::double precision)) AND _st_dwithin((g1.geom)::geography, (g2.geom)::g (...)"
"              ->  Seq Scan on parcels g1  (cost=0.00..20630.55 rows=152755 width=395)"
"              ->  Materialize  (cost=0.00..430.86 rows=7524 width=60)"
"                    ->  Seq Scan on powerplants_us_201603 g2  (cost=0.00..393.24 rows=7524 width=60)"


texts: Get the same coordinates for points after downloading a .osm file

texts:
I downloaded a .kml file that contains bike stations which I know the coordinates (latitude and longitude). So I got a layer that contains points corresponding to the stations. Then I downloaded the data to a  .osm file.

My problem is that the coordinates of the points (latitude and longitude) have changed in the .osm file, I can not distinguish the coordinates of  bike stations in this file knowing that the CRS used in QGIS is WGS 84; EPSG: 4326.

Why when downloading  the layer to .osm file QGIS changes the coordinates of the points?
What should I do to keep the same coordinates after downloading a .osm file please?


texts: Negative values in r.watershed flow accumulation

texts:
I am trying to process a DTM by using r.watershed in QGIS 2.14.3, with the GRASS GIS toolbox. It should be the file named: "Number of cells that drain through each cell".

I get some negative values, even if the shape of drainage line seems very correct. I think that the problem might be the number format (like float or integer).

I am working on a .tiff dem, but it was extracted by clipping from a .asc bigger dem.

How can I solve this?


texts:
I'm trying to cache OpenStreetMap related tile. The corresponding server has separate data related to opeenstreetmap (http://xx.xx.xx.xx/osm). It's an XZY Tiles service. Within the mapcache.xml file I set the source type to be tms service. But when I restart my Apache service I get a message:unknown source type tms for source osm.

Mapcache.xml file:

<?xml version="1.0" encoding="UTF-8"?>

<mapcache >
   <cache name="disk" type="disk">
      <base>C:/ms4w/tmp/ms_tmp/cache/primer3/{tileset}/{z}-{x}-{y}.{ext}</base>
      <!--<symlink_blank/>-->
   </cache>

   <source name="osm" type="tms">
      <getmap>
         <params>
            <FORMAT>image/png</FORMAT>
            <LAYERS></LAYERS>
         </params>
      </getmap>

      <http>
         <url>http://xx.xx.xx.xx/tile/{z}/{x}/{y}.png</url>
      </http>
   </source>

   <tileset name="test">
      <source>osm</source>
      <cache>disk</cache>
      <grid>g</grid>
      <format>PNG</format>
      <metatile>5 5</metatile>
      <metabuffer>10</metabuffer>
      <expires>3600</expires>
   </tileset>

   <default_format>JPEG</default_format>
   <service type="tms" enabled="true">
      <full_wms>assemble</full_wms>
      <resample_mode>bilinear</resample_mode>
      <format>PNG</format>
      <maxsize>4096</maxsize>
   </service>
   <service type="wms" enabled="true">
   </service>
   <service type="wmts" enabled="true"/>
   <service type="kml" enabled="true"/>
   <service type="gmaps" enabled="true"/>
   <service type="ve" enabled="true"/>
   <service type="demo" enabled="true"/>

   <errors>report</errors>
   <lock_dir>C:/ms4w/tmp/ms_tmp/cache/primer2</lock_dir>

</mapcache>

Why does it report that the unknown type is tms if it exists as a defined service?


texts: What does Error [sqlite3_rtree_geometry_callback is not found] mean in qgis 1.8 with windows 7 64 bit?

texts: qgis1.8 does not open. I get a message that sqlite3_rtree_geometry_callback is not found


texts: Creating in-memory raster layer using PyQGIS?

texts:
I have been researching quite a bit for creating in-memory raster layer using PyQGIS. i have done it through gdal "memory" layer. But i need QgsRasterLayer instance for displaying. I don't have the option of saving the results to disk and reading it again.
I tried doing this:

target_ds = gdal.GetDriverByName('MEM').Create('ras', x_count, y_count, 1, gdal.GDT_Int32)

target_ds= processing.runalg("gdalogr:rasterize",
                      {"INPUT": ellipse_shp,
                       "FIELD": "id1",
                       "DIMENSIONS": 0,
                       "WIDTH": x_count,
                       "HEIGHT": y_count,
                       "RAST_EXT": "%f,%f,%f,%f" % (x_min, x_max, y_min, y_max),
                       "TFW": 1,
                       "RTYPE": 4,
                       "NO_DATA": 1,
                       "COMPRESS": 0,
                       "JPEGCOMPRESSION": 1,
                       "ZLEVEL": 1,
                       "PREDICTOR": 1,
                       "TILED": False,
                       "BIGTIFF": 2,
                       "EXTRA": '',
                       "OUTPUT": "aa"})

But when i wanted to use it as:

target_ds["OUTPUT"]

it gave me an error as

  TypeError: str cannot be converted to qgis._core.QgsRasterLayer in this context

Although QGIS documentation says it supports all gdal sources, I couldn't make it work.

Can anyone help me?


texts:
I've excitedly just downloaded the latest version of QGIS 2.0.1-3 (QGIS Dufour) and now it won't open. I performed the exact same process last night (without problem) on my home laptop (64-bit OS & windows 7) using the correct installer package. Why wont this open now! Lisboa, (1.8.1 I think) worked fine. My computer is 64-bit OS and Windows 7.

Edits:

Have now tried the 32-bit version. Completely unresponsive. It went through an extra step during the install procedure loading a C++ package (from 2005). Other than that I'm still having the same problem.


texts: Get individual pixel values on a raster image using gdal and python

texts:
I'm now stuck with writing a piece of code to look up individual pixel values on a raster image.

This piece of code mainly executes two tasks:
(1) converting world coordinates to raster coordinates (I use a package named "Affine" to do that);
(2) then computing pixel offsets and figuring out which column and row the pixel is located at.

The major problem I'm facing now is that the code didn't return the correct column and row that the pixel is located at.

import sys, gdal, numpy
from affine import Affine

# coordinates to get pixel values for
lons = [-122.265373, -122.429139, -68.123095]
lats = [37.873090, 37.783640, 44.981935]

# open the raster file
ds = gdal.Open('data/nlcd_compressed.tif')
if ds is None:
    print 'Could not open the raster file'
    sys.exit(1)

# get raster size
cols = ds.RasterXSize #161190
rows = ds.RasterYSize #104424
bands = ds.RasterCount #1

# get georeference info
transform = ds.GetGeoTransform() # (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
xOrigin = transform[0] # -2493045.0
yOrigin = transform[3] # 3310005.0
pixelWidth = transform[1] # 30.0
pixelHeight = transform[5] # -30.0

# convert global coords to raster coords
aff = Affine.from_gdal(-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
x_coords, y_coords = aff*(numpy.array(lons), numpy.array(lats))
# converted coordinates
xValues = [-2496712.96119, -2496717.87417, -2495088.69285]
yValues = [3308868.8073, 3308871.4908, 3308655.54195]

band = ds.GetRasterBand(1) # 1-based index
#bandtype = gdal.GetDataTypeName(band.DataType) #Byte
#data = band.ReadAsArray(0, 0, cols, rows)
#value = data [122, 37] #127

data = band.ReadAsArray(xOffset)

# loop through the coordinates
for i in range(len(lons)):
    x = xValues[i]
    y = yValues[i]
    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)
    # get individual pixel values
    data = band.ReadAsArray(xOffset, yOffset, 1, 1)
    value = data[0, 0]
    print value


texts:
I badly need help...
I have a layer served from geoserver, this layer is in polygon.. and in every polygon the color varies according to a number. so it is like a choropleth map. now my problem is the color appears only at a certain zoom level in the map. When I zoom out the color assigned to the polygon disappears, then when I zoom in again the color comes back. What should I do to retain its color even in zoom in or zoom out.


texts:
I'm planning to create a map with GPS and routing support as well as other custom features for my local area. It will be for a commercial usage. I'm okay with showing OpenStreetMap logo in my app.

Is it okay to use OpenStreetMap data as base and add more features on my own for free?

Do I have to host my own server? (It will only support a map with around 4km square feet.)

But I do not wish to make what I created open-sourced.

Is OpenStreetMap for me?


texts:
I have been tasked to create a list of point symbols for the counties Watershed Improvement Program from scratch. They don't want modified icons from the programs database, but custom drawn. I already have the hand sketched symbols approved and I'm using Inkscape to create the vectors. Now to the actual question:
Is there a step-by-step process from converting the .bmp or .png files to Esri symbols? Ie. what size and format should the vector be created in, how do you upload into Arcmap (10.1), and how do you access your folder in the Style Manager? I know this may seem simple to a lot of you, but the Googled help I've seen is about as simple as Greek Calculus to me.


texts: Seeking GIS Library for C# Applications using WPF

texts:
This question has been asked multiple times before, but not in the last 5+ years so I think it's worth asking again as I'm struggling and the answer way well be different now.
I want to create a C# application (WPF) which has some GIS capabilities. I'm new to C#, application development and GIS so this should be an interesting challenge. My main objectives are:

Load (ideally free to use) OpenStreetMap tiles
Load Shape Files
Load DXF Files (not critical as could convert to SHP first)
Add points, lines and polygons
Interrogate SHP files
Allow a flat world projection for working in local coordinate systems.

For now the package will be used internally at work, but I wouldn't want to write off the possibility of selling it if it exceeds my expectations so I don't want to base it on free to develop but prohibitively expensive libraries if used commercially.
I have tried the following:

MapSui for WPF - Very easy to get a basic map started, but I find the documentation just stops after that (at least in a language I can understand) and you have to rely on their examples which for me are too involved for me to pick apart and re-create.

DotSpatial - I believe to be WinFrom only. Still I could use this.  Again, getting a basic map working was fairly simple (bit more involved than MapSui) but I couldn't work out how to only display the tools I wanted for example. Again, documentation seemed to stop after a basic set up.

ArcGIS Runtime WPF - Documentation appears to be excellent but it's not free. I find their Pricing a little confusing. It appears that if you don't use their locator services then it may be free? They mention a monthly developers charge, but unclear if this is for additional services that I don't need? The terms also say you have to allow location services and you can't use it for the "Runtime in a service" which makes no sense to me. If anyone can clarify the above that would be useful. If it can be used for free with local data that would be great.

SharpMap - Hasn't been updated since 2014 so I may avoid.

GMap - I don't believe it can operate in flat world or load SHP files which could make it's use very complicated.

ThinkGeo Desktop Maps - $1000 per year which isn't inline with my expectations from this package.

I seek any suggestions or thought on the best route forward for me?


texts: Calculating number and sum of length for linear objects inside polygon

texts:
I have two shapefiles.

the first includes Polygons. Each polygon has a unique "FID"
the second has linear objects that belong exclusively to one of these polygons and they also have a "CLASS" attribute which takes values 1 to 6 and a "LENGTH" attribute

For each of the polygons I need to calculate the number of linear objects that have a class value = 1, the number of linear objects that have a class value = 2 and so on.

Also for each of these classes that belong to each polygon I need the sum of their length.

For example, in the end i need to have in the polygon shapefile these values, i.e. each polygon will have 12 new values:

"FID", "Count1", "Count2", ..., "Count6", "SUM1", "SUM2", ..., "SUM6"

Is there any easy way to do this or do I have to try ArcPy?


texts: How to merge many shapefiles in different folders?

texts:
Directory structure is basically as such: county\City\prep\POI\Completed\C1, ...C2, ...etc\

I'm trying to write a script to go through the *\Completed*\POIS.shp and merge them all together, at a country level.  So the input would be the country folder as root.

So, to clarify:

county1\City1\prep\POI\Completed\C1\POIS.shp

county1\City1\prep\POI\Completed\C2\POIS.shp

county1\City2\prep\POI\Completed\C1\POIS.shp

etc

I started to write out the walk through the folders, but I figured I should get some guidance first.  So now I've got it working so if the merged files are in a specific folder I can easily run the XY coordinates and other selections I need to do.  The step just above where I've commented out should be the merge of all POIS.shp files and copied to the directory to run the other tools.  Hopefully (but probably not) not confusing.

import arcpy
import os
import sys

# Folder name to run XY script
folder = arcpy.GetParameterAsText(0)
##countryFolder = arcpy.GetParameterAsText(1)
# workspace
arcpy.env.workspace = folder

# List shapefiles
arcpy.AddWarning("Building List of Shapefiles...")
fcList = arcpy.ListFeatureClasses()

## Merge a Country's completed POI files (Country\*\PDA\Completed_Field_Files\*\POIS.shp\)
##for dirs in os.walk(countryFolder)
##    arcpy.Merge_management("POIS.shp",

# Add coordinates # code is good and works!
try:
    for fc in fcList:
        arcpy.MakeFeatureLayer_management(fc, "temp_" + fc)
        shapeFile = "temp_" + fc
        # Add XY coor to table
        arcpy.AddXY_management(shapeFile)
        arcpy.AddWarning("Building" + " " + shapeFile + " " + "XY Coordinates Successful")
        # Select only records we want with SQL: Action<>D; Position=Yes,NA; POI_ST_NUM<>blank
        arcpy.AddWarning("Selecting" + " " + shapeFile + " " + "Records...")
        try:
            arcpy.SelectLayerByAttribute_management(shapeFile,"NEW_SELECTION","POI_ST_NUM <> ' '")
        except:
            arcpy.AddWarning("No Blank Street Numbers")
        try:
            arcpy.SelectLayerByAttribute_management(shapeFile,"REMOVE_FROM_SELECTION","ACTION_COD = 'D' OR ACTION_COD = 'FR' OR ACTION_COD = 'AFR'")
        except:
            arcpy.AddWarning("No Action Codes to Remove")
        try:
            arcpy.SelectLayerByAttribute_management(shapeFile,"REMOVE_FROM_SELECTION", "POSITION = 'NO' OR POSITION = ' '")
        except:
            arcpy.AddWarning("No Position Errors to Remove")
        # Copy selected records in feature claseses to new shapefile
        arcpy.AddWarning("Creating New Shapefile...")
        arcpy.CopyFeatures_management(shapeFile, "RWT_XY_" + fc)
    arcpy.AddWarning("***********************")
    arcpy.AddWarning("Huzzah! File/s Complete")
    arcpy.AddWarning("***********************")
except:
    arcpy.AddError("*********************")
    arcpy.AddError("DOH! Adding XY Failed")
    arcpy.AddError("*********************")

Thanks!  By the way, this is my 2nd script ever with no previous training or knowledge, so be easy on me.  :-D

EDIT: try to make it more clear

I have dozens POIS.shp files I want to merge, but they are all located in different folders according to the structure I listed above. Always in a C1, C2, C3, ..., C16 folder.  The C folders are always in Completed_Field_Files folder.  The structure may differ a little bit, so I want the tool to crawl and merge all POIS.shp files within the Completed_Field_Fieles\C# folders.  There's other POIS.shp files as well, so I ONLY want the ones within those specific folders.  I'm basically just merging a bunch of shapefiles that are located in different but specific folders


texts: Long range Trilateration from 3 lat, long, range

texts:
I've been looking at the following solution to a similar problem
Trilateration using 3 latitude and longitude points, and 3 distances
However with my problem the 3 points are spread over a large distance and the distance is a distance-over-ground (ignoring altitude). So the curvature of the earth comes into play.

Example of one of the data points is...[lat,long,km]

51.505348978413,-0.11270052661132,75.639168
51.4845086249581,-3.16947466601561,173.809152
55.9568033803799,-3.20266968478392,465.100416

This can be visualised at http://www.freemaptools.com/radius-around-point.htm by pasting it into the section "CSV Upload - [latitude,longitude,radius(km)] per line"

Using the solution in Trilateration using 3 latitude and longitude points, and 3 distances I get the result..

51.8247121391 -0.802763718095

which isn't bad but it's about 14 miles out. I'm really wanting this down to ~5 miles of error.

I changed the ECEF calculation to one from here http://www.mathworks.co.uk/matlabcentral/fileexchange/7942-covert-lat-lon-alt-to-ecef-cartesian/content/lla2ecef.m and used an altitude of 0.

However, this is way way off the mark with...

56.2518456429 -0.813853131422

Intermediate results for xA,yA,zA...etc. are

lla2ecef

xA = [3978.79751502] yA = [-7.82628596545] zA = [5342.89431232]

xB = [3974.53105834] yB = [-220.086727248] zB = [5342.88794881]

xC = [3573.7981748] yC = [-199.973379186] zC = [5344.22411569]

original method

xA = [3965.56758195] yA = [-7.8002627162] zA = [4986.36680448]

xB = [3961.32002974] yB = [-219.355176278] zB = [4984.92406448]

xC = [3561.02859617] yC = [-199.258852046] zC = [5279.1109334]

So... the Z values for lla2ecef method are quite different. I'm guessing this is what's causing the issue, something to do with fixing the altitude to 0 and the fact my distances are distance-over-ground.

Am I going to need to do some projections instead?

My brain is about melted for today, can anyone shed light on ways to tackle this problem?


texts: How to resolve a rectangle with know points?

texts:
How do I find the coordinates(lat, lon) for D and E in the following rectangle?

I know:

A (lat, lon)
B (lat, lon)
C (lat, lon)
AE, BE, AC


texts: Projecting latitude and longitude such that the Euclidean distance is a good approximation of geographic distance

texts:
Goal:

I want to run various algorithms in Shapely on some geographic data that I have. It is acceptable to have some measure of error, as long as it is "reasonable". e.g. find closest pair.

Method:

I would like to project the points from latitude-longitude to cartesian space such that

euclidan_distance(P(p0), P(p1)) ~= geodesic_distance(p0, p1)

I was thinking of using PyProj but it may be overkill, and it doesn't seem the easiest thing to do.

Naive approach:

y = latitude * 110574
x = longitude * 11320 * cos(radians(latitude))

(numbers from https://en.wikipedia.org/wiki/Latitude#Length_of_a_degree_of_latitude )

Results of naive approach:

I put in a set of coordinates nearby, and got relative errors of 0.5%-15.3%. 15% seems excessive.


texts: Getting authentic and verified current location of user?

texts:
I am building an application that needs a user to send his current location. However, since GPS signals can be faked. I need the user current location to be authentic and verified.

I thought of using a Trusted third party like Google to verify and authenticate the user location before he can send it.

Is my solution based on a valid assumption?


texts: Geosoft GRD histograms do not save in layer file

texts:
My colleagues and I are using ArcGIS 10.3 and 10.1 with the free Geosoft plug in.
When saving layers containing geosoft grd files we usually display using a percent clip which require the histogram to be built. The issue is that despite saving a layer file the histograms need to be rebuilt for most (but not all) grids each time Arc is started.
Has anyone had this issue and do you know how to sort it out?


texts: Are there any Hosting companies, that offer Geoserver?

texts:
I am working with a non-profit, and we are planning to have a webmap for them. Most of the hosting providers just offer a simple LAMP setup on their Linux option, and if you want to install any software on your own, you need to go for a dedicated server (either real or virtual) which is very expensive for the organisation I am working with.

I have seen that there is a hosting provider for Mapserver based maps.

Is there any similar provider for Geoserver based maps?


texts: How do I create and manage custom SQL functions in CartoDB?

texts:
I'm working through a CartoDB tutorial at the moment, and got to the point of creating a custom SQL function.

Do I just write this in the "CUSTOM SQL QUERY" Field on the right of my table?

How do I edit previously created functions?


texts: Measure shortest distance from a point to a line in QGIS 3

texts:
I am fairly new to QGIS and I have been using it to determine the shortest distance between multiple points using a network of roads.  After seeing the results, I see that there is usually a gap between the origin point and the network (this happens with the destination points as well).  In some cases the shortest path cannot be determined is that gap is too large, but I do not know what is considered as "too large".
So my question is, is there a way to define a maximum allowable distance between a point and a line to determine a shortest path? (ie. If there is not a road within at least 200 meters from the origin then it is not possible to determine a shortest path between two points)


texts: Rectifying a Tiff Given One Corner Coordinate

texts:
I'm trying to rectify an image that I don't know much about other than the coordinate of the southwest corner, which I was able to pull using ExifTool.

Is there any programming solution to identifying the coordinates of the other 3 corners and using that data to warp the image?


texts: How to best create linked pdf to point data in ArcGIS server 10 sde 9.3?

texts:
I would like to create point data which offers a link-field to open a pdf file which I will have on my server. I have just read about "attatchments" which are new to ArcGIS 10, but this option isn´t available for my sde9.3 feature classes - only my filegeodatabases in 10.0 Desktop. I have sde 9.3 running with ArcGIS Server 10.0.

How would I best create pdf links to point features classes in SDE?


texts: Speeding up SQL spatial query gives spatial index error?

texts:
I am trying to split polylines by points which have a small buffer around them. Currently I have over 370,000 lines and 320,000 nodes and the query is running really slowly (I've left it for 3 days and it still hasn't completed). I have tried forcing a spatial index using with (Index(SI_tempPD)) but I get the following error:

  "The query processor could not produce a query plan for a query with a spatial index hint. Reason: Could not find required binary spatial method in a condition. Try removing the index hints or removing SET FORCEPLAN."

Below is the snippet of code that I'm trying to run when I get the error:

BEGIN INSERT INTO TempLines ( [linenum] ,[ogr_geometry] )
SELECT lines.[linenum] ,lines.[ogr_geometry].STDifference(points.[ogr_geometry].STBuffer(0.005))
AS ogr_geometry FROM dbo.TemplineData AS lines with(Index(SI_tempPD)) INNER JOIN dbo.[TemplineNodes]
 AS points ON lines.[ogr_geometry].STIntersection(points.[ogr_geometry]).STDistance(points.[ogr_geometry]) < 1
WHERE (lines.[linenum] <> points.[linenum]) END

Is there anyway i can speed up the query?

I also have a clustered primary key.  The execution plan shows that a filter takes up 36% of the cost and the insert takes up 64%

I'm using SQL Server 2008 (SQL server Management studio 10.50.1600.1)


texts: Encoding format in QGIS, Vietnamese accents

texts:
What's the data source encoding to use in QGIS to delete the Vietnamese accents?

I want to get my data in the attribute table from Bình Chánh to Bing Chanh.


texts: Implementing Pairwise Intersect Tool in ArcMap?

texts:
This is the first time I have tried to add a tool to the ArcMap 10.4.1 Toolbox.
I downloaded the Pairwise Intersect Tool from here and installed it in My Toolboxes in the ArcCatalog : http://www.arcgis.com/home/item.html?id=120807ac7f8449f89840ae48e87d7f55

I gave the following parameters, which I thought I understood from here:
http://pro.arcgis.com/en/pro-app/tool-reference/analysis/pairwise-intersect.htm

in_feature1 : Table
in_feature2 : Table
out_feature_class : Feature Class
join_attributes : String
Cluster_Tolerance : Linear unit
Output_type : String

But when I run the script I get this error

Executing: pairwiseintersect C:\Users\gree1117\Desktop\Tonya\Wytham\flight_path2_convex_hull.shp C:\Users\gree1117\Desktop\Tonya\Wytham\Wytham_landuse.shp C:\Users\gree1117\Desktop\Tonya\Wytham\fp2_convexhull_landuse5.shp ALL "5 Meters" INPUT
Start Time: Wed Mar 08 16:30:44 2017
Running script pairwiseintersect...
Wed Mar 08 16:30:45 2017
Processing features...
Failed script pairwiseintersect...

Traceback (most recent call last):
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcToolbox\Scripts\PairWiseIntersect.py", line 166, in <module>
    pairWiseIntersect(inputFC1, inputFC2, outputFC, fldList2Transfer)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcToolbox\Scripts\PairWiseIntersect.py", line 123, in pairWiseIntersect
    for row2 in cursor2:
RuntimeError: A column was specified that does not exist.

Failed to execute (pairwiseintersect).
Failed at Wed Mar 08 16:30:47 2017 (Elapsed Time: 3.01 seconds)


texts: Batch merging shapefiles from different folders, different subfolders

texts:
I tried some of the scripts I found but none worked so far.
What I am trying to do is batch merging hundreds, maybe thousands of shapefiles from different folders/sub-folders into a single shp or feature class (I read is better in case of large data).
My data is stored in [Shape Files] folder, as you may see here:

O:\ABC\DEF\2010\project_nu\Shape Files*.shp
O:\ABC\DEF\2011\project_nu\Shape Files*.shp
O:\ABC\DEF\2012\project_nu\Shape Files*.shp
O:\ABC\DEF\2013\project_nu\Shape Files*.shp
O:\ABC\DEF\2014\project_nu\Shape Files*.shp

(where there's many project_nu in each year)

I don't have any experience in Python just in ModelBuilder.

I did a little bit of research and I understand that first I need to use the os.walk through all these folders and search for *.shp polygons (all shp's are polygon so there's no need to test this). Then apply merge to the generated list.

I tried ModelBuilder too (Iterate feature class-> Collect Values-> Merge) and got stuck in the iterative process that kept running with no results (I even tried for a few shps).


texts: Performing selection on Attribute and Location of two different feature classes using ArcGIS for Desktop?

texts:
I have a point FC and a Polygon FC. Both FC have the attribute VILLAGE.
The point with Village = 1500 should be completely within the Polygon with Village = 1500.

How can I make a selection not only based on the location, but also on the attribute? In the End I want to see all Points/Polygons where the VILLAGE does not correspond.


texts:
I am working on an interactive map which I envision to have a similar structure to diffusion networks. I have data to pin point an event and I would like to show those involved with said event, and display this data year by year.

These networks are normally shown in abstract space and can be created in statistical software. I was wondering if there is any map services or current programs that allow for a geographic analysis of this sort.

Some examples:

http://thebritishgeographer.weebly.com/the-spread-of-disease-and-its-management.html
(under HIERARCHICAL DIFFUSION)

Diffusion Networks and Cascade Analytics
http://people.tuebingen.mpg.de/diffusion-networks-wsdm14/


texts: ENVI 5: High resolution export of multiband imagery

texts:
I have to export some multiband Landsat imagery (Landsat channel combination 543) in high resolution in ENVI 5 for a Figure which I later will edit in Adobe Photoshop and Illustrator, best would be TIFF (geoinformation is not important).

My attempts to do this:

I stacked all channels into one file with LayerStacking and then loaded the needed band combinations into a new view. However when I am trying to save that file, I can only select the stacked file with all bands.
I know that there is the Chip View to File option but the output
resolution is relatively low.

Is there a way to export it in higher resolution or another way?


texts: What is the meaning of Altitude in EXIF of an image?

texts:
Currently I have some photographs captured by DJI, and in their EXIF I saw that it contains: GPS Altitude, Absolute Altitude, Relative Altitude, GPS Latitude, GPS Longitude.

These images will be used for creating TIFF file by Pix4D. After creating TIFF file from those images, I use the GPS Lat/Long of one of those images and then extract the elevation from TIFF file and that elevation is different with the altitude from the image.  I also use Google Elevation API and it also different from the altitude of the image (and the elevation from my raster, too).

For example, with this lat/lon (GPS value from EXIF) 15.99849675/108.246482777778 I get those 3 values:

Altitude (GPS and Absolute) in the EXIF: 62.9m Above Sea Level
Relative Altitude: 49.80
Elevation I get from TIFF: 21.9400272369
Elevation I get from Google API: 4.45

So my questions are:

What is the meaning of the altitude in EXIF of an image?
The accuracy between the elevation extract from TIFF and the Google elevation API


texts: The position and size of my Tile Map is incorrect when it is displayed in React-Leaflet

texts:
I converted the GeoTIFF file to a tiles map. I'd like to display it in react-leaflet, but the size and position are not right.
What can I do to fix it?
This is condition when displayed.
(overlaid on OpenStreetMap whose position is correct)

If working properly, it should be located in about 100-meter square a place in Japan.
This is my TileLayer component for displaying my tiles map.
<TileLayer
    url="http://localhost:3333/assets/tms/{z}/{x}/{y}.png"
    noWrap="true"
    tms="true"
/>

Here is my gdalinfo from my GeoTiff file.
Driver: GTiff/GeoTIFF
Files: GeoTiff.tif
Size is 9105, 5977
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (138.638412499999987,35.633808000000002)
Pixel Size = (0.000000210840198,-0.000000163275891)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 138.6384125,  35.6338080) (138d38'18.28"E, 35d38' 1.71"N)
Lower Left  ( 138.6384125,  35.6328321) (138d38'18.28"E, 35d37'58.20"N)
Upper Right ( 138.6403322,  35.6338080) (138d38'25.20"E, 35d38' 1.71"N)
Lower Right ( 138.6403322,  35.6328321) (138d38'25.20"E, 35d37'58.20"N)
Center      ( 138.6393723,  35.6333201) (138d38'21.74"E, 35d37'59.95"N)
Band 1 Block=9105x1 Type=Byte, ColorInterp=Red
  NoData Value=-10000
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=9105x1 Type=Byte, ColorInterp=Green
  NoData Value=-10000
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=9105x1 Type=Byte, ColorInterp=Blue
  NoData Value=-10000
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=9105x1 Type=Byte, ColorInterp=Alpha
  NoData Value=-10000

Here is my tilemapresource.xml from my TileMapService.
<?xml version="1.0" encoding="utf-8"?>
    <TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0">
      <Title>GeoTiff.tif</Title>
      <Abstract></Abstract>
      <SRS>EPSG:4326</SRS>
      <BoundingBox minx="138.63841249999999" miny="35.63254780813838" maxx="138.64033219999999" maxy="35.63380800000000"/>
      <Origin x="138.63841249999999" y="35.63254780813838"/>
      <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
      <TileSets profile="raster">
        <TileSet href="0" units-per-pixel="0.00001349377265" order="0"/>
        <TileSet href="1" units-per-pixel="0.00000674688633" order="1"/>
        <TileSet href="2" units-per-pixel="0.00000337344316" order="2"/>
        <TileSet href="3" units-per-pixel="0.00000168672158" order="3"/>
        <TileSet href="4" units-per-pixel="0.00000084336079" order="4"/>
        <TileSet href="5" units-per-pixel="0.00000042168040" order="5"/>
        <TileSet href="6" units-per-pixel="0.00000021084020" order="6"/>
      </TileSets>
    </TileMap>

I tried too epsg:3857, the same CRS as OSM, but it didn't change.


texts: I have 7 shapefiles which I would like to combine using the 'Merge shapefiles' function in QGIS.  I got the error message that the shapefiles must all be the same geometry.  I checked this, one shapefile has the geometry 'WKB type - PolygonZ' - I think this was a shapefile created when I used the Vector-dissolve function.  All the other files have the geometry 'Polygon'. How do I change the geometry of this shapefile - or is it a case of recreating and re-saving the original shapefile?


texts: Choosing cell size in LAS Dataset to Raster Tool

texts:
I am creating a raster from an LAS Dataset (LAS Dataset represents only the ground return).

I am wondering how do you choose the sampling size? (which I have set to cell size)

My point spacing is 1.886 and my point density is 0.2811.