texts: Add data from a REGION FeatureCollection to my polygons FeatureCollection

texts:
I have generated a FeatureCollection of polygons and calculated the centroid and added it. I now want to pull in some data from FAO/GAUL/2015/level2 dataset with fields such as ADM2_CODE etc. This can be just the centroid.

How do I get data out of one FeatureCollection and add it to mine?

var fromList = ee.FeatureCollection("FAO/GAUL/2015/level2");

var filter = ee.Filter.inList('ADM2_NAME', ['Gwent']);
var filteredArea = fromList.filter(filter);

Map.setCenter(1.1, 51.5, 8);
Map.addLayer(filteredArea, {}, 'Regions');

var gfc2018= ee.Image('UMD/hansen/global_forest_change_2018_v1_6').clip(filteredArea);

var lossImage = gfc2018.select(['lossyear']).gte(7);
var lossBand = lossImage.and(gfc2018.select(['lossyear']).lte(20));
var zones = lossBand.updateMask(lossBand.neq(0));

// Add the loss layer in red.
Map.addLayer(lossImage.updateMask(lossBand),
            {palette: ['FF0000']}, 'Loss 2007 onwards');

var area_threshold = 300; // about 30 hectares?

// Define whole county as a region
var geom = filteredArea.geometry();

// Convert the zones of the forest loss areas to vectors.
var vectors = zones.addBands(zones).reduceToVectors({
  geometry: geom,
  crs: lossBand.projection(),
  scale: 30,
  geometryType: 'polygon',
  eightConnected: true,
  labelProperty: 'zone',
  reducer: ee.Reducer.count()
}).filterMetadata("count","greater_than", area_threshold);

// Make a display image for the vectors, add it to the map.
var display = ee.Image(0).updateMask(0).paint(vectors, '000000', 1);
Map.addLayer(display, {palette: '000000'}, 'vectors');

var SRTM = ee.Image("USGS/SRTMGL1_003"); // exampleelevation image

var slope = ee.Terrain.slope(SRTM);

// add min, max and mean elevation of each geometry
var reducers = ee.Reducer.median()
                .combine(ee.Reducer.min(), 'Elevation', true)
                .combine(ee.Reducer.max(), 'Elevation', true)
                .combine(ee.Reducer.mean(), 'Elevation', true);

vectors = SRTM.reduceRegions(vectors, reducers, 30);

// This function computes the feature's geometry area and adds it as a property.
var addCentroid = function(feature) {
  return feature.set({CentralPoint: feature.geometry().centroid(10)});
};

// Map the area getting function over the FeatureCollection.
var CPAdded = vectors.map(addCentroid);

// Get Admin Districts
var AreaName = filteredArea.aggregate_first("ADM2_CODE")

// How can I add the ADM2_CODE or any of the fields to my featurecollection? It can be just for the central point which I have calculated?

// Export the FeatureCollection to a CSV file.
Export.table.toDrive({
  collection: CPAdded,
  description: 'ForestPlus',
  fileFormat: 'CSV'
});


texts: No such file or directory PERMANENT\\MYNAME from Python+Grass?

texts:
I have an existing Python script that processes raster images with Gdal / Numpy. Now I want to extend the script with Grass. The python script is developed on Windows 7 with Eclipse/Pydev. The script is executed from within Eclipse.  I installed Grass with the OSGeo4W installer (Version 6.4.3).    I’m new to Grass, but as far as I understand I first need to create a new GRASS location. I try to do this with:

grass.core.create_location("C:\grassdata", “test”, epsg=21781)

however this throws the following error:

IOError: [Errno 2] No such file or directory:'C:\\grassdata\\test\\PERMANENT\\MYNAME'

create_location first runs g.gisenv, then g.proj and then again g.gisenv and finally tries to open C:\grassdata\test\PERMANENT\MYNAME that does not exists. Actually nothing was created in C:\grassdata.  This is the complete script:

import os, sys
import numpy as np
import gdal

def get_env(env):
  value = ""
  if env in os.environ:
    value = os.environ[env] + os.pathsep
  return value

input_data = r"C:\data\dom_eth\DOM_Clip.tif"

# start configure grass
gisbase = os.environ['GISBASE'] = r'C:\OSGeo4W64\apps\grass\grass-6.4.3'

os.environ['GISRC'] = r'C:\Users\Username\AppData\Roaming\GRASS6\grassrc6'
grass_dirs = [os.path.join(os.environ['GISBASE'], i) for i in ['etc', os.path.join('etc', 'python'), 'lib', 'bin', 'scripts']]
path = get_env('PATH') + os.pathsep.join(grass_dirs)
os.environ['PATH'] = path
os.environ['GRASS_SH'] = r'C:\OSGeo4W64\apps\msys\bin\sh.exe'

if not os.path.exists(gisbase):
  os.mkdir(gisbase)

sys.path.append(os.path.join(os.environ['GISBASE'], "etc", "python"))
import grass.script as grass
import grass.script.array as garray
import grass.script.setup as gsetup

gisdbase = r"C:\grassdata"
location = "test"

grass.core.create_location(gisdbase, location, epsg=21781)

A grep -rl "MYNAME" over the source code of Grass did not reveal something useful to me.
What is the correct way to use Grass with Python from the "outside"?


texts: Using output of iterate files as input for merge tool in ModelBuilder?

texts:
I am using trying to use an iterator to run through and folder of shapefiles and add them as input to the merge tool in ArcMap ModelBuilder. Something seems amiss; when I connect the iterate output file to the merge tool in ModelBuilder, it doesn't seem to take.

What am I doing wrong here?


texts: How to scatter plot the raster data of LST and NDVI?

texts: From the scatter plots, the dry edge and wet edge is to be determined to calculate soil moisture. I have used the MODIS products of LST and NDVI for determining the soil moisture index.


texts:
Windows Server 2016

ArcGIS Desktop 10.4.1

I've installed ArcGIS 10.4.1 on Windows Server 2016. When I connect to the server through Remote Desktop, ArcMap crashes on launch. ArcCatalog will launch while connected in RDP. However, I can launch ArcMap without any issues when physically on the server. I can open ArcMap on the server and then remotely connect and it will continue to function, but it will not launch a new instance.

I've tried deleting the normal template, updating the registry for the Python 2.7 installation, uninstalling Python, etc. None of those should be the problem since it works fine outside of RDP.

I don't know much about dump files, but the Exception Information is "The thread tried to read from or write to a virtual address for which it does not have the appropriate access." So, assuming it's trying to access something on the client machine? I've tried checking and unchecking all of the local resources.

I'm pretty sure I initially installed it while connected through RDP, but I uninstalled and reinstalled while not in RDP and that didn't change anything either. Not sure if I've missed something with Remote Desktop on the server.

Any ideas?

I've hit a dead end and didn't get an answer on the Esri forums.


texts: How put order in multiple raster file in RasterMosaic in FME?

texts:
I have couple of raster files which I want to mosaic in a specific order.
How I can do that?


texts: Geocode plentiful internationals adresses with Postgres or QGIS

texts:
I've got a database with a list of events (about 30 000 events) all around the world.
I would like to geocode those. For this, I have differents column (city/state/country).
I tried to geocode with QGIS plugin MmQGis, but I've got to much requests...
Then I wanted to try with tiger geocode, but it seems work just for US adresses.

Do you have any solution to geocode those internationals events ?


texts: ArcObjects: How to insert data into a table?

texts:
I have a table in ArcCatalog named WorkOrderAss.

This table contains 3 columns: (OBJECTID, FeatureName, Name).

I want to insert data to this table from the code behind in C#.

Any help would be appreciated.


texts: How to add X/Y data from Excel file (saved as CSV) on shapefile?

texts:
I need to add my X/Y coordinate data as points from my Excel I need help with adding points (using X/Y coordinates). Can anyone help?

import arcpy
from arcpy import env

env.workspace = "C:/myFolder"
env.overwriteOutput = True
out_path = "C:/myFolder"

out_name = "X_Y_Points.shp"
geometry_type = "MULTIPOINT"
has_m = "DISABLED"
has_z = "DISABLED"


texts:
Is there a way to directly add a WFS layer in ArcGIS Pro?

Prior to ArcGIS Pro 2.0 the only way to use a WFS service is to convert it to a feature layer.


texts:
I've installed R on Ubuntu 19.10 and successfully installed the plugin 'Processing R Provider'.  However, when I try running a script (e.g. Scatterplot) I get an error (see the report at the end of the post).  I have queried this with the developers of the plugin, but it would appear to be a local issue.

I can launch R from a terminal as well.

Any suggestions?

2019-09-29T16:32:51     WARNING    Traceback (most recent call last):
              File "/usr/share/qgis/python/plugins/processing/gui/ProcessingToolbox.py", line 239, in executeAlgorithm
              ok, message = alg.canExecute()
              File "/home/neil/.local/share/QGIS/QGIS3/profiles/default/python/plugins/processing_r/processing/algorithm.py", line 279, in canExecute
              msg = RUtils.check_r_is_installed()
              File "/home/neil/.local/share/QGIS/QGIS3/profiles/default/python/plugins/processing_r/processing/utils.py", line 315, in check_r_is_installed
              universal_newlines=True) as proc:
              File "/usr/lib/python3.7/subprocess.py", line 775, in __init__
              restore_signals, start_new_session)
              File "/usr/lib/python3.7/subprocess.py", line 1522, in _execute_child
              raise child_exception_type(errno_num, err_msg, err_filename)
             PermissionError: [Errno 13] Permission denied: '/usr/lib/R'

2019-09-29T16:38:57     WARNING    Traceback (most recent call last):
              File "/usr/share/qgis/python/plugins/processing/gui/ProcessingToolbox.py", line 239, in executeAlgorithm
              ok, message = alg.canExecute()
              File "/home/neil/.local/share/QGIS/QGIS3/profiles/default/python/plugins/processing_r/processing/algorithm.py", line 279, in canExecute
              msg = RUtils.check_r_is_installed()
              File "/home/neil/.local/share/QGIS/QGIS3/profiles/default/python/plugins/processing_r/processing/utils.py", line 315, in check_r_is_installed
              universal_newlines=True) as proc:
              File "/usr/lib/python3.7/subprocess.py", line 775, in __init__
              restore_signals, start_new_session)
              File "/usr/lib/python3.7/subprocess.py", line 1522, in _execute_child
              raise child_exception_type(errno_num, err_msg, err_filename)
             PermissionError: [Errno 13] Permission denied: '/usr/lib/R'


texts:
I am trying to create ellipses from fields in my posgres table. I found this custom Ellipse function and have been trying to tweak it to meet my needs.

Create Ellipse (x,y,rx,ry,rotation,#of segments in 1/4 of ellipse):

CREATE OR REPLACE FUNCTION Ellipse(double precision, double precision, double precision, double precision, double precision, integer)
RETURNS geometry AS
$$
SELECT ST_Translate( ST_Rotate( ST_Scale( ST_Buffer(ST_Point(0,0), 0.5, $6), $3, $4), $5), $1, $2)
$$
LANGUAGE 'sql';

My inputs into this are a lat and lon. Smajor/Sminor measurements in meters. I receive my rotation in degrees but have a written a function that has converted them to radians per the requirements as stated in some of the other forums I've seen with this similar problem. I have tweaked the above method to this:

CREATE OR REPLACE FUNCTION Ellipse2(double precision, double precision, double precision, double precision, double precision, integer)
RETURNS geometry AS
$$
SELECT ST_Translate( ST_SetSRID( ST_Rotate( ST_Scale( ST_Buffer( ST_SetSRID( ST_Point(0,0), 3395), 0.5, $6), $3, $4), $5), 4326), $1, $2)
$$
LANGUAGE 'sql';

I create the point (0,0) and then set it's SRID to 3395(world meters based system) in hopes that the follow on functions use meters when it extrapolates the smajor/sminor sides. I'm getting ellipses that are massive though (i.e. continent size not 600 meter size) I'm pretty sure my issue is with projection based units. Can someone point out to me what I am doing wrong?


texts:
For one of my projects using PostGIS 2 I have a table with GPS tracks (igc_files). It holds records including a MULTIPOINT() column that has about 1,000 points for each record and a column that has an array of the corresponding timestamps. I also have a raster table (elevations) that should at some point hold all of the SRTM elevation data. For now I have only imported a few relevant areas though.

What I want to do is get the elevation for each point of one certain record and the corresponding timestamps. I have implemented an naive approach and it does work, but I'm wondering about its performance implications. This is how it currently looks:

SELECT timestamps[location_id] AS timestamp,
       ST_VALUE(elevations.rast, subq.location) AS elevation
FROM elevations, (
    SELECT (ST_DUMP(igc_files.locations)).path[1] AS location_id,
           (ST_DUMP(igc_files.locations)).geom AS location
    FROM igc_files
    WHERE igc_files.id = 1
) AS subq, igc_files
WHERE igc_files.id = 1 AND
      ST_INTERSECTS(subq.location, elevations.rast) AND
      ST_VALUE(elevations.rast, subq.location) IS NOT NULL

Since this is a GPS track and the points are not that far apart, PostGIS could assume that the next point is in the same raster tile and check that tile first before checking the index again from the top. Is that actually done by PostGIS or is it just starting to search the index from the top for each of the points?

I have tried to use ST_INTERSECTION() instead of ST_VALUE() and feeding in the MULTIPOINT directly, but that was actually even slower by a few orders of magnitude. Is there a feature/function of PostGIS that I am missing to do this more efficiently?


texts: Copy and paste values from selected features to a table in ArcMap?

texts:
I am trying to find a more efficient way to select features from a polygon feature class and copy over their GUIDs to a separate table in ArcMap.  I have a layer of pavement polygons and I would like to select all the sections of pavement that are receiving maintenance treatment, and then copy over their values into another separate table.  The separate table has a 1-Many relationship with the pavement polygon layer.

Is there a tool in ArcMap or a scripting option to achieve this?

Here is a visual representation of my problem:

These are the selected polygon features and their attribute table.  I want to copy the GUID values over to this table:

Note how this allows there to be a 1 Pavement Polygon to Many treatment records in the table possible.  Right now I just do them one by one, but I was hoping for another solution where I could copy/paste over many polygon GUIDS into the table at the same time.


texts: Superposing two plots with Geopandas, with non fill colours in one of them

texts:
I would like to plot a map of the edges of the French departments, and the heat maps at the lower scale of the French IRIS. I read both shapefiles using GeoPandas.

And now

dep.plot(color="white", edgecolor='black')

and

iris.plot(column = 'ratio',cmap = 'Purples',ax=ax)

I would like to visualise the black edge colours of the departments on the IRIS plot, but when I try this

fig, ax = plt.subplots(figsize = (20,16))
df.plot(column = 'ratio',cmap = 'Purples',ax=ax)
dep.plot(color=None,edgecolor='k',linewidth = 4,ax=ax)

the first figure does not appear, and is hid by the second. I have tried to use alpha, to put the colour as None, but nothing seems to work.

Any help?


texts: "Regular points" and "Vector grid" don't work for me in QGIS 2.18.4

texts:
I have a QGIS 2.18.4. I wanted to use "regular points" feature, but whenever I try QGIS gets stuck at the beginning

Nothing happens afterward.

The same thing happens with another research tool: "vector grid".

I have no issues with other QGIS components, except this two.

Is there anything I can try, except of uninstalling QGIS and trying another version?


texts: What are the ways we can control the quality of a DTM?

texts: I imagine that for DTM obtained from Stereoscopic procedures or interpolation from a set of topographic points a comparaison with the original data can help, but what for DTM obtained from Remote Sensing ?


texts:
The picture is an example, but the method would be the same at mine bigger dataset.
So i have some parcels with each a ID value and some fields.
such as field 1, field 2, field 3, fieldd 4,... with some numbers in the columns.

Now the objective is to calculate the average for each parcel (ID), so for each row and locate the answer in a new field 'average'.

how is it possible? i only know the way to calculate a average for a column, not for a row...


texts: How to fill a hole in a multi-polygon layer on QGIS 2?

texts:
I need to fill a hole in a multi-polygon layer on QGIS 2.0.
I've set the snapping option for the layer as follows :

mode: to vertex and segment
tolerance 4 pixels
avoid intersection "enabled"
topological editing "enabled"

When I try to edit a rough shape around the hole and confirm it right clicking, the system told me "The feature cannot be added because its geometry collapsed due to intersection avoidance"

Update:

I have a multipolygon file with a hole. I'm not able to fill the hole because if I manually edit it the system tells me that it is impossible to add it because the geometry is collapsed due to intersection avoidance. If I try to use the specifc tool it tells me that it's impossible due to the snapping settings.

These are the images:


texts: Selecting multiple features with rectangle using ArcGIS API for JavaScript?

texts:
I am using ArcGIS API for JavaScript 3.7.

I am trying to see if I can select features in the same fashion that shift+click and dragging lets you zoom into the map by drawing a rectangle. The trouble that I am encountering is that I want this functionality on a certain layer of the map, but I do not know if this layer is ready yet. I've tried using the following:

map.on('load', function initSelectToolbar(event) {
    selectionToolbar = new Draw(event.map);
    var selectQuery = new Query();

    on(selectionToolbar, "DrawEnd", function (geometry) {
        selectionToolbar.deactivate();
        selectQuery.geometry = geometry;
        featureLayer.selectFeatures(selectQuery, FeatureLayer.SELECTION_NEW);
    });
});

I would like to see if it is possible to run this code when I know for sure that the layer is ready. As far as I am aware, I can't change the feature layer the code references once it is loaded in.

EDIT:
It looks like I got it to work. I found a way to grab the layers URL using map.getLayer(layerId).url (how I missed that is beyond me), created a new feature layer for selection, and removed it when the user was done. I still think this could be done better though.


texts: Seeking Symbology for Interactive Mapping?

texts:
When a user clicks on a map what should they expect to see?

I would say a symbol for where they clicked.  A pushpin seems to be the most accepted icon for generic data. Google has pushed the upside down tear drop.

What are other ideas/options and is there a place for finding professional gis oriented icon packs?

I ended up creating this,

This just popped up from the guys at map box. Clean maki icons


texts: Adding SQL layer with dropdown selector in Carto?

texts:
I am trying to make a map that has a layer showing District boundaries of South Africa, with a selector allowing the user to view which climate change indicators are a priority to each District. Using Carto's latest docs (with the Builder), I have successfully managed to build the foundation of my map and add the first layer of the District boundaries (see screenshot below):

Now I am trying to use the dropdown menu (id="selector") to search through a second dataset (called indicators) and fill the polygons that answer "yes" to the relevant feature column. Below is the database I created using QGIS:

Here is my code for the SQL layer:

    //Add SQL layer
    const indicatorsSource = new carto.source.SQL(`
      SELECT *
        FROM indicators WHERE selector ilike 'yes'
    `);
    const indicatorsStyle = new carto.style.CartoCSS(`
      #layer {
        polygon-fill: #162945;
        polygon-opacity: 0.5;
      }
    `);
    const indicators = new carto.layer.Layer(indicatorsSource, indicatorsStyle);

    //Add layers to map
    client.addLayers([boundaries, indicators]);

I think using "selector" in the SQL query is where I am going wrong.

The full code with a working SQL example can be viewed here


texts: Identifying/exporting adjacent polygons in ArcGIS Desktop?

texts:
I have many polygons (US Counties) with agricultural land data in them. I would like to identify and export (ideally into excel) the data from the counties that are adjacent to the ones with agricultural land use data.

Is there a way to do this quickly or en masse in ArcGIS Desktop 10.2?


texts: Extract values from raster, R and ArcGIS different results

texts:
My goal is to extract values from the raster called pftnc by using the other raster sp after converting it to a point shapefile. Files available here pftnc and sp

I have two raster files

pftnc=raster("pftnc.tif")
sp=raster("sp.tif")
sp
    class       : RasterLayer
    dimensions  : 142, 360, 51120  (nrow, ncol, ncell)
    resolution  : 96486, 96515  (x, y)
    extent      : -17367529, 17367529, -6356742, 7348382  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : Puma_concolor
values      : 0, 1  (min, max)

pftnc
    class       : RasterLayer
    dimensions  : 720, 1440, 1036800  (nrow, ncol, ncell)
    resolution  : 0.25, 0.25  (x, y)
    extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
    value     : 0,1

In R, first I use aggregate to upscale one to have the same dimensions, I crop and then I project.

pftagg <- aggregate(pftnc, fact=4)
pftext <- c(xmin= -180, xmax= 180, ymin= -60, ymax= 90)
pftagg=crop(pftagg,pftext)
pft1 <- projectRaster(pftagg, sp)

Then I convert raster to points and extract the values

sppoints=rasterToPoints(sp,fun=function(x){x==1})
result=extract(pft1,sppoints[,1:2])

I do the equivalent operations in ArcGIS 10.5, but I don't reproject or crop because when I import the layers everything lines up perfectly.

So the problem now is that I have inconsistencies between the two methods. Aside from a few missing points from the ArcGIS operation, I get slightly different means and standard deviation. I also tried to import in ArcGIS the raster after cropping, aggregating, etc, and I still get these differences. So things in ArcGIS are at least consistent with itself. I am wondering if these small differences are due to how the two software are dealing with projecting the latlong raster to CEA projection.


texts:
I have a query to calculate how many crimes (of a certain type) are rural and urban:

select st_intersects(c.wkb_geometry,u.wkb_geometry) as ur,count(c.*)
from crime.street_crime c, meridian.urban as u
where c.crime_type = 'Vehicle crime'
group by ur;

This takes several seconds to run on a smallish dataset, explain says:

HashAggregate  (cost=147319.21..147319.73 rows=2 width=7043)
  Group Key: ((c.wkb_geometry && u.wkb_geometry) AND _st_intersects(c.wkb_geometry, u.wkb_geometry))
  ->  Nested Loop  (cost=74.84..144643.18 rows=535205 width=7043)
        ->  Bitmap Heap Scan on street_crime c  (cost=74.84..2772.55 rows=2893 width=233)
              Recheck Cond: ((crime_type)::text = 'Vehicle crime'::text)
              ->  Bitmap Index Scan on street_crime_type_idx  (cost=0.00..74.11 rows=2893 width=0)
                    Index Cond: ((crime_type)::text = 'Vehicle crime'::text)
        ->  Materialize  (cost=0.00..41.77 rows=185 width=6810)"
              ->  Seq Scan on urban u  (cost=0.00..40.85 rows=185 width=6810)

so it is ignoring my spatial index, that I expected it to make use of.

If I rewrite the query to:

select st_intersects(c.wkb_geometry,u.wkb_geometry) as ur,count(c.*)
from crime.street_crime c,
(select st_union(wkb_geometry) as wkb_geometry from meridian.urban) as u
where c.crime_type = 'Vehicle crime'
group by ur;

it takes miliseconds, explain says:

HashAggregate  (cost=3587.76..3588.28 rows=2 width=265)
  Group Key: ((c.wkb_geometry && (st_union(urban.wkb_geometry))) AND _st_intersects(c.wkb_geometry, (st_union(urban.wkb_geometry))))
  ->  Nested Loop  (cost=116.15..3573.29 rows=2893 width=265)
        ->  Aggregate  (cost=41.31..41.32 rows=1 width=6810)
              ->  Seq Scan on urban  (cost=0.00..40.85 rows=185 width=6810)
        ->  Bitmap Heap Scan on street_crime c  (cost=74.84..2772.55 rows=2893 width=233)
              Recheck Cond: ((crime_type)::text = 'Vehicle crime'::text)
              ->  Bitmap Index Scan on street_crime_type_idx  (cost=0.00..74.11 rows=2893 width=0)
                    Index Cond: ((crime_type)::text = 'Vehicle crime'::text)

This makes no sense to me, as by unioning the geometry I am effectively defeating any spatial indexing as the majority of my points will fall in the bounding box and need a more expensive point in polygon check but it is much faster.

Why doesn't my spatial index make the naive query work fast?


texts: when using arcpy.mapping.ListDataFrames, python.exe crash

texts:
Has anybody had this problem? this drive me crazy. what is wrong with my code?

import arcpy
mxd = arcpy.mapping.MapDocument(r"D:\Project\Project.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]

python crash at the last line


texts: Remove unconnected polylines from a shapefile ArcPy

texts:
How do you isolate a group of polylines that are not connected to the bigger polyline group.

Ideally, I would like to rank all the polyline groups based on how many features there are in each group.

The arrow in red shows the smaller polyline group that has to be isolated.

I have used,

arcpy.SelectLayerByLocation_management(temp, "BOUNDARY_TOUCHES", temp, "", "NEW_SELECTION", "INVERT")

But, this only works if the "group" has only one polyline.

I could also use the code below but it takes hours to run

arcpy.SelectLayerByAttribute_management("temp", "NEW_SELECTION", "FID = 100")
arcpy.SelectLayerByLocation_management("temp", "BOUNDARY_TOUCHES", temp, "", "ADD_TO_SELECTION")
arcpy.GetCount_management("temp")


texts: Loop over pixel data using Python in ArcMap

texts:
Further to my question on SO I'm trying to check the validity of bitmaps. Mainly for transparent or black pixels.

I'm new to ArcMap, and I now know that Python scripts can be run from Arc.

Any hoops, the following script loops over the bitmap pixel by pixel using Pillow.

# get image (theFile) & threshold (th)
f = theFile
img = Image.open(theFile)

w = img.size[0]
h = img.size[1]

pix = img.load()

blackCount = 0
for y in range(h):
  for x in range(w):
    p = pix[x,y]
    r,g,b, a = p

    print ("r, g, b: %s, %s, %s" % (r, g, b))

    if (r <= th) & (g <= th) & (b <= th):
      blackCount +=1

However, I realised that Pillow can't handle large bitmap files (I'm looking at images in excess of 2GB)

What's the best way to achieve this through ArcMap?

Bearing in mind I need to look over hundreds of directories containing each tiff. The results just need to be in the form of a text file - damage to each bitmap(s) (number of pixels effected) and the name of the file(s) in question.


texts:
EDIT: I've tried running the code linked in the duplicate flag, it's posted here for reference, and it doesn't run. I had already referred to that answer in my initial post.

My problem

Downloaded shapefile is missing .shx file and trying to use pyshp to build it on a mac.*

I've tried this code that I've picked up from the pyshp documentation. When checking the r.shapes and r.records, it looks like all of the info is there.

myshp = open("Stockton/Stockton.shp", "rb")
mydbf = open("Stockton/Stockton.dbf", "rb")
r = shapefile.Reader(shp=myshp, dbf=mydbf, shx=None)

In the documentation, it looks like each record and shape is being added one at a time. Is there a method I'm missing? Or could a list comprehension be used here to loop through everything?**

Data and References

I've downloaded the Microsoft Building Data for Stockton CA from here. I've tried to follow the code from the answer here

    # Build a new shx index file
    #Code by Joel Lawhead http://geospatialpython.com/2011/11/generating-shapefile-shx-files.html
    import shapefile
    # Explicitly name the shp and dbf file objects
    # so pyshp ignores the missing/corrupt shx
    myshp = open("Stockton/Stockton.shp", "rb")
    mydbf = open("Stockton/Stockton.dbf", "rb")
    r = shapefile.Reader(shp=myshp, shx=None, dbf=mydbf)
    w = shapefile.Writer(r.shapeType)
    # Copy everything from reader object to writer object
    w._shapes = r.shapes()
    w.records = r.records()
    w.fields = list(r.fields)
    # saving will generate the shx
    w.save("myshape")

but get this error:

TypeError                                 Traceback (most recent call last)
<ipython-input-94-345df1656b96> in <module>
      7 mydbf = open("Stockton/Stockton.dbf", "rb")
      8 r = shapefile.Reader(shp=myshp, shx=None, dbf=mydbf)
----> 9 w = shapefile.Writer(r.shapeType)
     10 # Copy everything from reader object to writer object
     11 w._shapes = r.shapes()

...

TypeError: expected str, bytes or os.PathLike object, not int

*Because I'm using a Mac, none of the solutions here seem like they would work.

**Still fairly new to Python so please forgive wrong/confusing terminology.


texts:
I'm trying to upgrade multiPrint plug in to QGIS3, but I'm newbie on python

My upgrade return an error like:

AttributeError: 'QgisInterface' object has no attribute 'activeComposers'

activeComposer no more exist in QGIS3

The code:

def populateComposerList(self):
    #called to populate field list for WHERE statement
    wdgt=self.dlg.composerList
    wdgt.clear()
    for cView in self.iface.activeComposers ():
        item=QListWidgetItem()
        item.setFlags(item.flags() | Qt.ItemIsUserCheckable)
        item.setCheckState(Qt.Unchecked)
        item.setText(cView.composerWindow().windowTitle())
        wdgt.addItem(item)
    wdgt.sortItems()


texts: Save output of Con statement to different folder

texts:
I have a folder of shapefiles that I want to convert to rasters and then apply a Con statement to them. I set my workspace to where the shapefiles are and do Point to Raster conversion and Con with the following code:

import arcpy
import os
from arcpy import env
from arcpy.sa import *

#Set workspace to Max Temps folder and allow for overwrite
env.overwriteOutput = True
env.workspace = r"C:\Users\uma.bhandaram\Desktop\2015\DailyMax2015"
outputFolder = r"C:\Users\uma.bhandaram\Desktop\2015\RasterCon2015"

#Get Max Temp shapefiles and convert them to Raster
shpFileList = arcpy.ListFiles("*.shp")
for shpFile in shpFileList:
    shpFileName = os.path.splitext (shpFile)[0]
    rastFile = shpFileName + ".tif"
    arcpy.PointToRaster_conversion(shpFile, "Max_2015", rastFile, "MOST_FREQUENT", "", 0.05)

#Get raster files and run through Raster Calculator to see if Max threshold has been met
rastFileList = arcpy.ListFiles("*.tif")
for rastFile in rastFileList:
    rastFileName = os.path.splitext (rastFile) [0]
    MaxRstIn = Raster(rastFile)
    MaxRstOut = Con(MaxRstIn, 1, 0, "VALUE <=85")
    MaxRstOut.save("Con" + rastFile)

How do I save the output rasters of the Con function in a separate folder? I have specified an output folder but am stuck on pointing my Con files to go here.


texts: Stray horizontal lines after converting from GeoJSON to topojson

texts:
After converting a GeoJSON file to topojson, and loading the topojson into Leaflet, stray lines appear in the resulting map.

Like this:

What is causing this? If I view the GeoJSON, these stray lines do not appear. I've played with the quantization behavior, and this does not seem to help. thanks.

(I am using Leaflet-Omnivore to digest the topojson.)


texts:
I did all that the qgis page says. Install ELGIS. after that i investigate and tried different commands.. but nothing works..the last I tried was:

yum --enablerepo=elgis install qgis

this was the response:

Loaded plugins: fastestmirror, refresh-packagekit, security
Loading mirror speeds from cached hostfile
 * base: mirror.edatel.net.co
 * epel: mirror.symnds.com
 * extras: mirror.edatel.net.co
 * rpmfusion-free-updates: mirror.hiwaay.net
 * updates: mirror.edatel.net.co
Setting up Install Process
No package qgis available.

seems just with centos this is complicated!


texts: Installing SDA4PP plugin for Kriging in QGIS

texts:
I am trying to do kriging in QGIS. I have loaded the SDA4PP plug in.
Almost all of the items in the SDA4PP dropdown menu are not highlighted and are unresponsive when clicked on - including the kriging option.

Any ideas?


texts: KML 2.2 parsing with a group of point Geotools

texts:
I am doing a viewer with geotools and I am trying to read a kml file. I am using PullParser class for parsing and also I am using KML and KMLConfiguration classes from org.geotools.kml.v22 package. If I use a kml file with just one geometry, the code works, but when I try with a file with more than one point, it read only the first coordinates, even though it capture all the names of the points.
I don't know how to make possible to read more than one point of the same file
Do you know how can I solve this?

Here is the code:

void readKML(File file){
 fis = new FileInputStream(file);
 PullParser parser = new PullParser(new KMLConfiguration(), fis, KML.Placemark);
 int index = 0;
 DefaultGeographicCRS crs = DefaultGeographicCRS.WGS84;
 SimpleFeatureTypeBuilder typeBuilder = new SimpleFeatureTypeBuilder();
 typeBuilder.setName(file.getName());
 typeBuilder.setCRS(crs);
 typeBuilder.add("name", String.class);
 typeBuilder.add("the_geom", Geometry.class);
 SimpleFeatureType type = typeBuilder.buildFeatureType();
 SimpleFeatureBuilder builder = new SimpleFeatureBuilder(type);
 DefaultFeatureCollection newCollection = new DefaultFeatureCollection();
 StyleMap.GeomType geom = null;
 SimpleFeature f = null;
 while ((f = ((SimpleFeature) parser.parse())) != null) {
   String geometryDesc = f.getDefaultGeometry().toString();
   //This function is for making the style
   geom = getGeometryType(geometryDesc);
   Geometry geometry = (Geometry) f.getDefaultGeometry();
   Object nameAttribute = null;
   nameAttribute = f.getAttribute("name");
   builder.addAll(new Object[]{nameAttribute, geometry});
   SimpleFeature feature = builder.buildFeature(String.valueOf(index++));
   newCollection.add(feature);
  }
}

I based on this reference:
How to parse KML data using geotools?
and also on this code:
http://direct.massapi.com/source/github/17/47/1747678994/plugins/org.locationtech.udig.catalog.kml/src/main/java/org/locationtech/udig/catalog/kml/core/KmlUtils.java.html#85


texts: Select Polygon Intersections Between Two Tables With No Common Key

texts:
I am new to geodatabase queries and cartodb.
I am trying to write a SELECT query that will do the following.

I have two tables.

Table 1 contains congressional district polygons and info about representatives.

Table 2 contains service territory polygons and info about providers.

There is no common key between the tables.

For each polygon in Table 1, I need to find all of the intersecting polygons from Table 2, and return all of the calculated polygon intersections.

Help?


texts: Returning gps points in form of gpx files in pgrouting?

texts:
I am trying to create a service that returns the shortest path for a given start and end destination.  What I did was used shortest path query provided by pgrouting to calculate the vertex and edges that form the shortest path. Then for each of the edges I calculated their corresponding coordinates and send them as points in a gpx file. Now when I plotted the coordinates and joined them by a straight line, it gave me a path but not the curve that I needed.

I am wondering what is the best way to send that information, I mean the path information. Currently I am just sending the coordinates. So lets say for a long curve path, do I need to actually send all the intermediate points in the curve to actually plot it like a curve or is there other way.

What is the standard procedure?


texts: Using ArcGIS to filter cells where 7 or 8 of neighbors are different

texts:
I have land cover data that I would like to clean using some generalization tools. Majority filter does almost exactly what I want, but is missing one key option. Majority filter allows you to change cells where 5 out of 8 of the neighbors are a different class. The problem is that in my data there are many instances when 5 out of 8 is too aggressive. I would like to only target those cells that are surrounded by 7 or 8 different neighbors.

Is there a way to do this in ArcGIS?


texts: Geotiff files - Zoom into a specific zipcode/latlong/apn

texts:
I am new to Geo tools.
I am trying to assess the DSM/DEM of buildings in specific cities. So far, I could obtain a geotiff file from https://earthexplorer.usgs.gov/ and plot maps of the entire earth using rasterio in python following this tutorial: https://www.earthdatascience.org/courses/use-data-open-source-python/intro-raster-data-python/fundamentals-raster-data/
After much research, I am unable to find an easy way to zoom to the granularity of a single apn/house and plot its elevation.


texts:
I use a method based on the "project a point" sample to project points from S-JTSK (wkid 102067) to Web Mercator (wkid 102100) and they are ~100m away from their correct position. This happens when I convert them to WGS84 (wkid 4326) as well. Points in other coordinate systems are projected correctly. I tried to use another geometry server, but the problem is obviously not there. Explicitly setting the transformation didn't change anything and setting forward transformation is not suitable for projected to projected SRS. When transforming to WGS84, I end in the Atlantic and the point even doesn't draw there. I tried to use WGS84 as an intermediate step, but the final result is the same.

My code:

reprojectStart: function (point)    {
    var sr = new esri.SpatialReference({
         wkid: thisDijit.map.spatialReference.wkid
    });
    var gsvc = new GeometryService("https://tasks.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer");

    var params = new ProjectParameters();
    params.geometries = [ point ];
    params.outSR = new SpatialReference(102100);
    params.transformForward = true;
    var def = gsvc.project(params, function(projectedPoints) {
        var pt = projectedPoints[0];
        point.x = pt.x;
        point.y = pt.y;
        point.spatialReference = pt.spatialReference;
    });
    def.then(function(results){
        return this.anotherMethod(point);
    });
}

The only workaround I know now is to either avoid S-JTSK or have both map and the points in it.

I wonder whether it is a bug ("coordinate system not fully supported") and I should report it to ESRI, or if I have some error in my code. Btw. it's a Web AppBuilder custom widget, so it's possible that it works otherwise but clashes with something in WAB.

EDIT: Transformation wkids I've tried: 1623,1625, 4827, 4829 and 5239. @mkennedy's comment encouraged me to try them again and I found the transformation is not being set correctly. params.transformation.wkid causes params.transformation not defined error. params.transformation = {wkid:5239} makes an object, but AFAIK true transformation object should have more attributes than just wkid. Other ProjectParameters attributes are correct. So the question basically changes to how to set the transform parameter when the way in the official example doesn't work?

EDIT2: I started examinating another way, through proj4js. However, I face some problems while making it work.


texts: What is the smallest size a point can be represented on ArcMap 10.2?

texts:
I have a large number of points that I need to plot and want to avoid overlaying them wherever possible. Whenever I choose a figure <1 for their size, I don't see any difference, yet ArcMap seems to accept the figure.

Is size '1' the smallest they can be represented?


texts:
I'm trying to find an R package to duplicate a specific bit of functionality of the ESRI network analyst: What ArcMap refers to as 'Service Areas'. Related terms include

catchment areas
isochrones
'reachability' (walkalytics package)
network buffers

What I want looks like an isochrone, but rather than measuring the amount of distance covered in a specified unit of time, I want a polygon representing the bounds of all points within 1.25 miles over the street network.

Nothing in the open streetmap routing machine (osrm package) meets my
needs.
Walkalytics isochrone would give me what i need, if I knew what
the walking speed was.

Could someone either point me to a package with this functionality, or tell me what the Walkalytics default speed is?

I'd really like to be able to do this in R, rather than having to go back into ArcMap everytime I need to determine the network distance reachable from a given metro station.


texts:
I have a large line feature class (15GB) in a File Geodatabase that I would like to clip.
I am trying to use ogr2ogr to do this, however I keep getting the error:

ERROR 1: Unable to write feature 11802 from layer inputLines.
ERROR 1: Terminating translation prematurely after failed
translation of layer inputLines (use -skipfailures to skip errors)

If I use –skipfailures, a lot of my lines are missing from the output.

The feature class itself seems ok.

The ArcToolbox's Check Geometery tool returns no errors with the
input feature class.
I am also able to use ogr2ogr to copy the input
feature class without any errors and have tried to use new feature class as input with the same result.

Here’s what I’m using when I try to output to another feature class.

ogr2ogr e:\output.gdb e:\input.gdb  inputLines -clipsrc clipPoly.gdb -nln outputFc -f "FileGDB"

Here’s what I use when I try to output to PostGIS

ogr2ogr PG:"dbname='mydb' host='localhost' port='5432' user='postgres' password='password'" e:\input.gdb  inputLines -clipsrc clipPoly.gdb -f "PostgreSQL"

The error is nearly identical except for the feature id number.

Last, I've also tried the ArcToolbox Clip tool, it just crashes after about 10 minutes.


texts: Snap points shapefile to line shapefile using shapely

texts:
I have 2 datasets: 1) a road network shapefile which has many polylines and 2) a bus stop shapefile which has many points representing bus stops. My points shapefile has CRS of EPSG 4326 and my line shapefile has CRS of EPSG 29617. I am not sure if I should and how to convert them.

I would like to use Python only and available packages (e.g. shapely, fiona, ogr) to snap the closest bus stop that is within 45 metres to the nearest road. Please note and see picture: there are multiple bus stops (green dots) near a road(brown line). Only the closest green dot (i.e. one green dot only) must be snapped to the brown line.

I have written this code based on many forums which have parts of the answer. I have never coded in Python before.

# path of inputs
line = "C:\folder\line.shp"
point = "C:\folder\point.shp"

# open line and points data
gdf_segments = geopandas.read_file(line)
gdf_hec_points = geopandas.read_file(point)

#join all line and points
shply_line = gdf_segments.geometry.unary_union
shply_point = gdf_hec_points.geometry.unary_union

# perform interpolation and project point to line
pt_interpolate = shply_line.interpolate(shply_line.project(shply_point))

The result of the code should be a new point shapefile where the points have been moved to the line.

Reference: I found the tool in ArcGIS (http://pro.arcgis.com/en/pro-app/tool-reference/editing/snap.htm) but I cannot use any software like QGIS, GRASS, ArcGIS etc. The code needs to be in Python only.


texts:
I am fairly new to ArcGIS Pro.  I have created an offline Mobile Map using ArcGIS Pro 2.3.0. which I am developing for a waypoint mapping project in Android.  The Mobile Map uses the Imagery basemap. The map localised to a small area in Caloundra, QLD, Australia.

The issue I have is when I enter a waypoint used Lat. Long. coordinates it appears on the Imagery map offset about 3 meters south-west.  I wish to shift the map north-east to align the view of where the waypoint marker (a Graphic on a GraphicsOverlay) should be (in relation to the imagery).

I have tried a few approaching including creating a map using 'GDA94 / Geoscience Australia Lambert (3112)' spatial reference.  However, I have yet to resolve this issue.

So my question is, how to I adjust (shift) a map?


texts: Connect to single table within a MSSQL database or terminate connection to database once layer added to QGIS project?

texts:
I'm connecting to MSSQL server '123' that contains 100s of databases and I am connecting to just database 'ABC'. 'ABC' contains roughly ~700 tables.  Once I add the 'xyzTable' as a layer to my QGIS project, the connection to the ~700 other tables in 'ABC' remains.  This is putting unnecessary resource strain on the 'ABC' database and '123' server.
Is there anyway to only connect to the table(s) I need from 'ABC'? Or is there a way to terminate the connection to all the other ~700 tables in 'ABC' without closing out QGIS and reopening the project.


texts: Handling raster data in a "rotated pole projection" with R?

texts:
My question is about using R to transform raster data from a rotated pole projection to a regular latlon grid - and then plotting it.

I am working with data from a numerical weather model, which is specified in a rotated latlon grid - I do this entirely in R. The data is provided in grib format, which is easily read with "raster" package. As example data, I have a grib file in a format identical to the weather data but with ground elevation values for my area of interest (Denmark and southern Sweden) - (example data available here).

Loading and plotting the data looks like:

library("raster")

geopotential <- raster(dk_geopotential_file)

plot(geopotential,
     xlab="Longtitude (rotated)",ylab="Lattitude (rotated)")

This looks good - shapes of Denmark and Sweden is clearly seen. Asking for the crs information (as a proj4 string) of the rotated grid gives:

> crs(geopotential)
CRS arguments:
 +proj=longlat +a=6367470 +b=6367470 +no_defs

Now I want to transform this into a regular WGS84 lat-lon format of the type:

"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

I would normally do this with the projectRaster function:

geopotential_wgs84 <- projectRaster(geopotential,
                                    crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")

But nothing happens! The data remains in the rotated latlon grid.
I have tried following these two threads (here and here), and played around with "ob_tran" arguments (Oblique Transformation).

Choosing CRS / PROJ4 string for EURO-Cordex rotated pole projection? seems to do this for a shapefile with the sf package, but I can't make it work with my raster data.

In my case, the coordinates of the geographic South Pole (not the North Pole as in the two linked threads) in the rotated grid is: (lon;lat) = (3;-40).


texts: Map rotation field don´t change in Atlas QGIS 2.18.14. Is this a bug?

texts:
I set up Scale and Rotation Map Fields to use with Data defined Override function, in an Atlas.

As I browse the maps in the Atlas, the Scale field appears fine, it change with the data from that field, but the Map rotation field don´t change, and has data in it.

Is this a bug?

I try with several types of fields, and still don´t show the data in it...


texts: Custom C# tool populate field combobox parameter from layer parameter

texts:
I'm creating a custom tool for arcmap toolbox where user select FeatureLayer in one combobox, and in next combobox user can select a field from selected FeatureLayer.
FeatureLayer combobox fills automaticaly from:

...
inputParameter.DataType = (IGPDataType)new GPFeatureLayerTypeClass();
inputParameter.Value = (IGPValue)new GPFeatureLayer();

The problem arises when I need to fill field combobox. I'm not realy sure how to set up its parameters. right now i'm trying something like this:

inputParameter = new GPParameterClass();
inputParameter.Enabled = true;
inputParameter.Name = "FieldTrg";
inputParameter.DisplayName = "Feature Class Fields";
inputParameter.Direction = esriGPParameterDirection.esriGPParameterDirectionInput;
inputParameter.ParameterType = esriGPParameterType.esriGPParameterTypeRequired;
inputParameter.DataType = (IGPDataType)new GPStringType();
inputParameter.Value = (IGPValue)new GPString();
//inputParameter.AddDependency("FeatureClassTrg");
//IGPFeatureSchema fgs = new GPFeatureSchemaClass();
//fgs.FieldsRule = esriGPSchemaFieldsType.esriGPSchemaFieldsAll;
//inputParameter.Schema = (IGPSchema)fgs

In UpdateParameters function I were able to get fields in selected FeatureLayer/FeatureClass, but I had no idea where and how should I put them back in combobox.
Also currently the field 'combobox' isn't really a combobox, but just a string field, so probably GPString isn't corect thing to use as its DataType/Value.


texts: Spatial analysis and parallel processing in R

texts:
I am trying to do a point-in-poly analysis using R. I have 4.5 million points and 5k polygons. This is computationally intensive. I was pointed to parallel processing through another post and tried to run point.in.poly through parallel processing.

Without parallel

table1 = point.in.poly(spdf, all_polygons)

With parallel

library(foreach)
library(doParallel)

registerDoParallel(makeCluster(2))
ptm <- proc.time()
foreach(i=1:length(split_df), .combine = rbind, .packages='spatialEco','sp')  %dopar%  {
  myvars = c("LONG_X","LAT_Y")
  producers_lat_long = split_df[i][myvars]
  spdf <- SpatialPointsDataFrame(coords = producers_lat_long, data = split_df,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
  clip1 = point.in.poly(spdf, all_polygons)
  getValues(clip1)
}
proc.time() - ptm
endCluster()

Without parallel results in about 7 minute processing time for a 100k subset. With parallel takes longer, then freezes and crashed my computer. Note that I am on a four core computer.

How might I speed up this function? Without parallel processing, my rough estimation is that this process could take 6 hours.


texts: Connecting directly to an Oracle Database using PyQGIS

texts:
I have an Oracle database that I would like to do INSERTS and UPDATES on.

I would like to do this using python inside a standard QGIS installation without having to install additional packages or clients.

Is this possible and in that case, how can it be done?


texts: Make Feature Layer won't display in Table of Contents

texts:
I have a point layer "SurfaceLocation".  It has a data field called "AgeRank".  One of the values it can be is "M".

In Modelbuilder, I want to take that layer, select only the points that have a value of "M", and create a new layer in the Table of Contents.  I'm using the "Make Feature Layer" tool.  I use "SurfaceLocation" as the input, apply the correct expression, name the Output layer "Minors" and set the output to "Add to Display".

I run this and nothing happens.  No new layer appears in the Table of Contents.  When I try running it again, I get an error saying the Output Layer "Minors" already exists.  It's not in the Table of Contents, so where is it?

I've tried closing and relaunching ArcMap.  I've tried rebuilding this in a new model (and it still says it already exists).  "Add results of geoprocessing operations to the display" is checked in geoprocessing options.

I've reviewed ModelBuilder not adding output feature class to Table of Contents? and decided to post this question because that didn't have an answer that worked for me.  I was hoping for answers that I could use.

I've tried making the output layer a parameter.  When I run the model, it doesn't even show up at the beginning for a value (the start parameter does).  Unfortunately, I don't know a thing about python yet.  I was hoping there was a tool that I couldn't find that would do this.

Further explanation:
Not sure why this simple problem can't have a simple solution.

The model is an oval (layer input) a rectangle (Make Feature Layer) and another oval (layer output).  I want the input to be an existing layer without any definition query.  I want the tool to apply the definition query and output another layer into the TOC.  I do not want to create another feature class. The model doesn't have to look like this.  If I need to string a couple of tools together to do it, that's fine.

* If I was doing this manually, I'd right click on the first layer, select copy, paste it, rename it, and add a definition query.  Then I would "Apply symbology from layer". *

I can't find a tool that does this.  "Make Feature Layer" doesn't put it in the TOC.  It puts it in_memory.  That would be fine if I could then find a tool that takes it out of memory and adds it to the TOC.  Being in_memory doesn't help me.  I want to display it, not do any processing to it.


texts:
I am very new to QGIS and have run into a problem dealing with shapefiles.
Recently, I went into the field with an app for smartphones called 'SW Maps' that I used to create some basic layers for lines and points. Under the 'point layer' I have three different types of points (I know I should have separated them when I created the layers).
Now  that I'm home, in front of the computer, I exported the shapefiles and loaded the data into QGIS.
Everything looks good, except for that I want to change the symbols and colors of the points so I can categorize them according to what they are.
The problem is that when I change the attributes to one point, it changes all the rest. I am not sure how to split(?) each individual point so I can edit them accordingly without having to create a duplicate, editing the color and symbol, making sure all the attributes copied over, and then deleting unwanted duplicates.
Does anyone know the best way to work around this? Or a good tutorial video on this type of editing? Or will I have to duplicate and edit individually?


texts: Outline overview for available datasets in QGIS

texts:
Assuming you have several 100 000 files on your discs

which are all georeferenced
which are all readable within QGIS (Raster maps, sqlite maps, vector maps, shapefile etc.)

After 20+ years in working with geodata I lost the overview over my collection.

In consumer tools like Oziexplorer or Oruxmaps you have a feature where you can see the outline of available datasets at your region of interest. If you click on this outline, the data is loaded.

Is there a tool like this in QGIS, maybe a plugin, and able to handle 100 000++ of datasets?


texts: What was the projection used for the Soviet Topographic Maps? http://maps.vlasenko.net/smtm50/k-37-072-1.jpg i am trying to georeference a map from Georgia. I believe the projection used was GK (3 degrees) with a Pulkovo datum (epsg 2530 in the case of the map above of the Georgian city of Poti)


texts: Geonetwork Empty Geographic bounding box

texts:
I set up Geonetwork 3.02 while creating a metadata geographic bounding box and spatial extent is empty. How can I fix these problems ? I use turkish language.
Can you help me please?


texts: Geoserver: how to call a Geotools function from Openlayers

texts:
I intend to create a custom Geotools function that does some special insertions into a PostGIS DB.

My question is how to call that function from Openlayers?


texts: How to count points once with overlappying polygons?

texts:
All,

We're using parcel data (polygon) to count student locations (point) based on parcel type. We're doing this using SQL and the ST_Intersects function.

It works great, except the City and County of Denver stores parcels with condominiums on them as multiple overlapping polygons, that is, each condo at a location gets its own parcel polygon - if there are 20 buildings on a parcel, there are now 20 coincident parcel features stacked on top of one another.

This is a problem because each parcel would count 1 student at that location once for every condominium parcel. So while there is only 1 student at that 'spatial location', if there were 20 condos at that location, that would be counted as 20 students.

Using a SQL 'distinct' does not seem to work...

How can we either:

a) Generalize the polygon data such that all overlapping polygons become a single polygon layer

b) Use some other SQL wizardry to eliminate the overlapping polygons from the query?

Thanks!!!

-m


texts: Shortcut to activate and deactivate a layer filter

texts:
I'm looking for the fastest way to filter and unfilter a layer in qgis 3.

My current workflow is too tedious: right click on layer / filter / write filter or delete filter

My final goal is to see at one sight the results of a filter/unfilter.


texts: Visualising a Raster (KAP open sea map file) using python with rasterio how do I get the colours correct

texts:
I have just started using rasterio to show raster maps in Python. I want to display some OpenSeaMap (KAP raster format). I have started with a single kap file, I can load and display a single kap file. However, the colours are not what I want. My code is below.
`
import rasterio
from rasterio.plot import show
src = rasterio.open("L10-296-560-8-24_10.kap",cmap='')
show(src)

`
This displays

But what I want is

Does anyone know how I tweek the colormap to get to what I want? Or maybe I should be using a different python library.


texts: OpenLayers Only draw features in a vector layer within view

texts:
I am dealing with a kml that I am importing into openlayers as a vector layer. This kml has many features so it dramatically slows down the web app. I would like to make it so only the features within a certain extent are drawn (based around the user's position and also only at a higher zoom level). So far I set up an event listener to check for the zoom level. Once the appropriate zoom level is reached the vector layer is drawn however, it draws the whole layer. Is there any way to do this such that it only shows part of the KML? I can't find the answer or an example anywhere.

map.events.register("zoomend", this, function(e) {
    var zoom = map.getZoom();
    if (zoom > 10){

           kml = new OpenLayers.Layer.Vector("Deer Accident Locations", {visibility: true,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
            url: "deer.kml",
            format: new OpenLayers.Format.KML({
                extractStyles: false,
                extractAttributes: true,
                maxDepth: 0

            }),
            style: {externalGraphic: 'images/Deer_green.png'  }          })
    });

    map.addLayers([kml]);


texts: Invalid topology error when running script?

texts:
I am trying to run a python script and I get an invalid topology rule when I get to the feature to line tool. I am using arc 10.1 but I think the script was written in 10.0. I know it was built in model builder and then converted. Invalid Topology [Feature not found on delete.] is the exact error.

Has anyone seen this before?

There are no topologies in the .gdb I am trying to process.


texts: Batch Update Fields in Field Calculator or Model

texts:
I have a dataset of block-level census data that I clipped to an Illinois State Senate district. To update the data (e.g. for population) I'm applying an area multiplier, i.e taking the area of each attribute after and dividing it by the area before and applying the percentage to each field. Since I wasn't aware of how to do a geoprocessing clip that would also clip data I'm left editing the table manually

My question is: how do I batch multiply columns against one column? I have 100 columns of data I need to multiply against a percentage contained in an aptly named proportion column. Can this be done through a model in ArcGis, or through a Python script in the Field Calculator? I have both QGIS and a demo license of ArcGis available.


texts: Problem with CSS Extension in Geoserver 2.7.0

texts: I am trying to load the CSS Extension in GS 2.7.0 by unzipping the contents of the archive into the WEB-INF/lib directory in GeoServer. When restarting GS i get a error : "Error creating bean with name "cssDemo" and a 503 Error while open the localhost


texts: OpenStreetMap coloring the map tiles by state, zip code polygons from C# VS 2013 WPF desktop application

texts:
My question is likely to close on stackoverflow.

My aim is to create a map from C# VS2013 WPF desktop application (not web application) with OpenStreetMap on win7.

This a colored map with US states colorized differently.

This can be shown in the attachment below.

I am new to WPF and handle the map tiles from XAML.

It seems to be very hard to find state or zip code borders from the map tiles and do coloring on them.

I have tried to add an image to the map to show colors but it cannot be aligned well with the adm areas/polygons. It was also stretched when zoom-in/out. I have tried to overlay one maptile on top of another one but there are always build errors. I have tried to search a similar  example code to find out how to implement the coloring for adm areas on map tiles, but I cannot find one.

The C# code in which I need to put an image (defined in "MyTileSource") on top of map tiles. I hope that the overlayed map tiles can be aligned well and zoomed/panned well. I got error:

  TileLayerCollection was not found in MapControl

But, actually, it has been defined.

This is only the first step. Currently, the "MyTileSource" is just a simple single-colored png. I still need to find out how to find state/county/zip-code borders in each map tile so that I can do coloring based on these administrative polygon areas.

 <map:Map >
  <map:TileLayers>
        <map:Map.TileLayerCollection>
            <map:TileLayer
                    SourceName="OpenStreetMap"
                    TileSource="http://{c}.tile.openstreetmap.org/{z}/{x}/{y}.png"
                    MaxZoomLevel="19"/>
            <map:TileLayer>
                <map:TileLayer.TileSource>
                    <map:MyTileSource/>
                </map:TileLayer.TileSource>
            </map:TileLayer>
        </map:Map.TileLayerCollection>
    </map:TileLayers>
</map:Map>


texts: Maintaining exact geometries with st_subdivide of PostGIS

texts:
I have noticed that when I use st_subdivide() function in PostGIS, the resulting area of the polygons can be different than the input.
Is there anything I can do to avoid this happening?


texts:
Does anyone know if there is some kind of documentation on which database types that you can connect to in QGIS? I would also like to know if:

Can I use the databes for read only or read & write.

thanks


texts: Overlapping polygons and determining true gaps

texts:
Using QGIS, I’m trying to merge several shp files into a larger area, the issue I’m finding is that some of these shp files contain “gaps” basically smaller polygons inside the larger one. I want to be able to bring in say 15-20 of these and want to have all the “gaps” that are overlapped by another polygon to be ignored and the ones that aren’t overlapped to remain. Here are some pictures of what I mean.

I don’t have trouble merging the polygons and cutting the gaps out of them to create a polygon with “holes” however I’m not sure how to go about telling the program to disregard “gaps” overlapped by other polygons.

If I’m only dealing with two or three polygons this is easily done by hand but I’m looking at doing 20 or more way more complicated polygons to create one larger one with only actual gaps left. I'm not sure if this can be done in the program or needs a script (I'm not good at writing scripts at all).


texts: pgAdmin does not integrate shp2pgsql in menu

texts:
Recently I updated my pgAdmin to 4, and I still want to use PostGIS Shapefile Import/Export Manager. I did not find any method to make this happen and even if shp2pgsql-gui still works outside pgadmin, this does not appear on pgAdmin's menu.

Am I doing something wrong, or is it no available on this version?


texts:
I have a GeoServer WMS with multiple layers, each with a time dimension (so its a WMS-T). I am able to query multiple layers simultaneously, but I am not able to query multiple TIME fields simultaneously. Is such a query possible, or would I have to make a script to query each TIME individually and then pool the results?

If this is possible with something other than Geoserver, or through some plugin, please let me know too.

This is the query I am trying to execute:

/geoserver/workspace/wms?SERVICE=WMS&\
VERSION=1.1.1&\
REQUEST=GetFeatureInfo&\
QUERY_LAYERS=workspace%3AACCGUST,workspace%3AT2M&\
LAYERS=workspace%3AACCGUST,workspace%3AT2M&\
INFO_FORMAT=application/json&\
FEATURE_COUNT=50&\
X=50&\
Y=50&\
SRS=EPSG%3A4326&\
WIDTH=101&\
HEIGHT=101&\
BBOX=-101.832275390625%2C31.42364501953125%2C-84.078369140625%2C49.17755126953125&\
TIME=2019-11-16T04:00:00.000Z,2019-11-17T04:00:00.000Z

And the output:

{"type":"FeatureCollection","features":[{"type":"Feature","id":"","geometry":null,"properties":{"Accum_Max_Gust":27.708324432373047}},{"type":"Feature","id":"","geometry":null,"properties":{"Temperature":33.7214469909668}}],"totalFeatures":"unknown","numberReturned":2,"timeStamp":"2020-01-07T18:09:04.642Z","crs":null}

As far as I understand I can't use a WFS because I am querying a large NetCDF file and the data is too big (when extracted) to store in a database.


texts: Making dropdown in ArcMap tool where each item links to spreadsheet?

texts:
The script below imports a spreadsheet into ArcMap and generates a shapefile. Currently, I can use the ArcMap tool that interfaces with the script to navigate to the spreadsheet I want to import. This tool will be used periodically to import from several spreadsheets, and I'd like to have the tool interface be just a dropdown, where each item links to a corresponding spreadsheet. Is this possible?

# Author: John K. Tran
# Contact: jtran20@masonlive.gmu.edu

import arcpy
import os
from os import listdir
import csv
print "arcpy imported"

"""
Input CSV can look like:

Lat,Lon,First Name,Last Name
12.34,56.78,Joe,Smith
0.98,7.65,Jane,Doe
65.43,43.21,Bob,Sagat

Make sure 'Lat' and 'Lon' fields are the 1st and 2nd column in the CSV
respectively.

"""

incsv = arcpy.GetParameterAsText(0)
#Change this to the path of your CSV file.

outfc = r'N:\ReferenceMaterials\ArcGIS Python\Pavement Cores by County\A
Interim Shapefiles' + '\\' + arcpy.GetParameterAsText(1) + '.shp'
#Change this to the path of your output FC.

spatialref = arcpy.SpatialReference(4326)
#Create the spatial reference object as WGS84. Can modify if desired.

testData = r'N:\ReferenceMaterials\ArcGIS Python\Pavement Cores by County\A
Interim Shapefiles'
if len(os.listdir(testData)) == 0:
    pass
else:
    for t in listdir(testData):
        if t.endswith('.shp'):
            arcpy.Delete_management(os.path.join(testData, t))

arcpy.CreateFeatureclass_management(os.path.dirname(outfc),
os.path.basename(outfc), "POINT", None, None, None, spatialref)

csv.register_dialect("xls", delimiter=",", lineterminator="\n")
#Register the dialect for native CSV syntax in Microsoft Excel.
f = open(incsv, "r")
reader = csv.reader(f, dialect = "xls")

headers = reader.next() #Read the first line as the header names.

for header in headers[2:]:
#Add fields for remaining columns if needed. Default is TEXT field.
    arcpy.AddField_management(outfc, header, "TEXT")

cursor = arcpy.da.InsertCursor(outfc, ['SHAPE@XY'] + headers[2:])
#Create InsertCursor.
count = 0
for row in reader:
    if count % 1000 == 0:
        print "processing row {0}".format(count)
    Ycoord = row[0] # Make sure 'Lat' is in the 1st column.
    Xcoord = row[1] # Make sure 'Lon' is in the 2nd column.
    newrow = [(float(Xcoord), float(Ycoord))] + row[2:]
    cursor.insertRow(newrow) # Insert point in FC for each row in CSV.
    count += 1

del cursor
f.close()


texts: Extracting pixel values by points from multi-band image and converting to table in Google Earth Engine?

texts:
I would like to extend this answer to Extracting pixel values by points and converting to table in Google Earth Engine? in order to extract values from a multi-band image.

I am using the example linked to manually select the band and export each GLCM texture measure as a 2D table (w/ values from each scene ID in columns).  I was hoping to run it as a single command (based on original answer from Rodrigo Principe).  The solution can either write a file for each texture (in a loop) or export a single, really wide table.

The collection is Sentinel-1 and the IW_H_OG is a feature collection of images over a study region (with the data converted from dB back to Original scale for GLCM calculations).

I've tried placing the  mapfunc  function in a for loop, but can't seem to get it right. I've also tried lopping the  .get("VH_diss");  off the end of the value variable (this might be working but GEE takes ages to export even a small test subset- so unsure if this is the simple solution).

I've pasted the tail end of the code below.  please excuse all the extra comments, I'm still coming to grips with JavaScript.

    //subset out only the GLCMs I want
        function subset(IMG){
          return IMG.select([0,1,2,3,4,5,8,14,18,19,20,21,22,23,26,32]);
          }
        var glcm_sub = IW_H_OG_glcm.map(subset);

    // create function that extracts pixel values for each point
    // from each image and writes to a table
    var mapfunc = function(feat) {
      var geom = feat.geometry();
      var addProp = function(IMG, f) {
        var newf = ee.Feature(f);

        // MY (BAD) WAY OF GETTING DATE INFO- STILL DEBUGGING
        var date = IMG.get('system:index');

        // SELECT BAND HERE- haven't figured out how to grab all
        var value = IMG.reduceRegion(ee.Reducer.first(), geom, 10).get("VH_diss");
        return ee.Feature(ee.Algorithms.If(value,
                                           newf.set(date, ee.String(value)),
                                           newf.set(date, ee.String('No data'))));
      };
      var newfeat = ee.Feature(glcm_sub.iterate(addProp, feat));
      return newfeat;
    }; // end of function

    // now apply function (with map) to our FeatureCollection (the cuts)
    // like ratio, these will be calculated and left on the OG scale
    var MY_pixels = myPoints.map(mapfunc);
    //print(MY_pixels.first());
    //01 VV_asm
    //02 VV_contrast
    //03 VV_corr
    //04 VV_var
    //05 VV_idm (homogeneity- Inverse Difference Moment)
    //06 VV_savg
    //07 VV_ent
    //08 VV_diss
    //09 VH_asm
    //10 VH_contrast
    //11 VH_corr
    //12 VH_var
    //13 VH_idm
    //14 VH_savg
    //15 VH_ent
    //16 VH_diss

    // export to drive
    Export.table.toDrive(MY_pixels,
    "Jam_III_2015_GLCM_16-16", //job name
    "test_GEE_export", //folder name
    "Jam_III_2015_GLCM_16-16"); //file name


texts: I need to create a circle using OGR and C# bindings.

texts: I have a reference point in latitude and longitiude. I need to create a circle in C# using GDAL and see if my list of points/polygons lie within the newly created circle. Is there a simple approach to doing this?


texts: Accessing ArcGIS server based services outside of firewall

texts: I would like to know the options for making ArcGIS based map services (ArcGIS Server 10) available to users who are out of our firewall, mostly field crews and users with web access in the field. These services are not secured at the moment. I am looking for some feedback from users who have implemented this. Thanks Jay


texts: Extracting values of one point dataset to second point dataset based on proximity using R or ArcGIS Desktop?

texts:
I have to point shapefiles and would like to extract values from one point shapefile to another one based on some measure of proximity, i.e. for any given point in dataset x, extract the mean value for all points in dataset y within a 10km radius of the given point in dataset x.  Or, for any given point in x, extract the value of the closest point in y and add to x.

I'm familiar with extract function in r to extract raster values to point. Not sure if sample data would help in this one, but would like to do this in r, or failing that, ArcGIS Desktop.


texts: Raster equivalent of "Recalculate Feature Class Extent"

texts:
Using Arc Desktop 10.7, I create an output raster using the Extract by Attributes tool.  However, the output raster retains the same x/y spatial extent as the input raster, even though its pixels occupy a different extent.

How do I recalculate the output raster's new spatial extent?


texts: Counting points in polygon for subset of points layer

texts:
I have a points layer with a column that separates the points into different types (e.g. 'high', 'medium', and 'low').
I need to count the number of points within polygons, for each of the different types.
At the moment, I'm doing it by creating a separate points layer for each of the types and then running the points within the polygon tool for each of them.
I was just wondering if there's a tool that lets you run points within polygon but with conditions so that I could run it for those points where the type is 'high' etc?


texts:
I am use to posting in Stack Overflow, but it seems this question is best geared towards the GIS Stack Exchange. I need  to understand the .ntf file types for satellite images. I am in grad school and studying about applications of ML, both supervised and unsupervised classification, to different interests. One of those areas of interest is satellite images. I have tried to do my best due-dilligence and research everything I can take in regarding images from remote sensors, including rasters/vectors, diffrent file types and processing techniques. A few of my references are:
USGS
NOAA
esri
ArcGIS

However, many of the papers I have read on ML and satellite image processing utilize .tif or .jpeg formats. I am having a very difficult time processing these .ntf files. From what I understand from Wikipedia , GDAL, & Purdue is that .ntf file-types are used by the government and are basically a grouping of multiple segments under one file:

from L3 Harris

Security Segments
Image Segments
Graphic/Symbol Segments
Label Segments
Annotation Segments
Text Segments Data
Extension Segments
Display Levels

I am hoping to train a simple unsupervised method (kmeans) to some of the images of this file type OR simply convert them, both of which are giving me trouble. For example, when using python I can utilize gdal to get a lot of good information, such as metadata, bands, and raster. However, all the images seem to have only one band and this is giving issues when trying to applying sklearn.cluster as the input is expected to be more than 1D.

Is there a different way to get the image data to be used in classification tasks?

Very simple for classification, just want to group by pixel - which, I am guessing is the issue since band 1 == black&white. I may be completely ignorant as well so someone do tell me to shut it.

I will include the current method below and perhaps someone can direct me to good resource for more enlightenment?

from osgeo import gdal, gdal_array
import numpy as np
from matplotlib import pyplot as plt

raster_img = gdal.Open('image_A.ntf')
band_n = raster_img.RasterCount
#band_n
#output > 1
raster01 = raster_img.GetRasterBand(1)
raster01_Arr = raster01.ReadAsArray()

#Trying to reshape for KMeans fails

classes = {'land':0,'water':1}
n_classes = len(classes)

rows, cols, bands = raster01_Arr.shape
##Error! 3 vars declared, only 2 values

rows, cols = raster01_Arr.shape

rows
#>4000
cols
#>4700

X = raster01_Arr.reshape(rows,cols,1) #added one for band

kM = KMeans(n_cluster = n_classes, random_state = 2).fit(X)
#fails, expects 2D, only received 1D

Example .ntf format files can be found here


texts: Iterate feature attributes of a postGIS layer using pyqgis

texts:
I am creating a QGIS plugin. Loaded postGIS layers onto QGIS. now, need to iterate the feature attributes with next and previous buttons.Once the layer is loaded, first attribute and feature is selected.

feature = selectedLayer.getFeatures().next()
 selectedLayer.setSelectedFeatures([feature.id()])

Using this code, first attribute and feature is selected. I am struck in selecting the next attribute of feature using next button.


texts: Disappearing Point Features in Geoserver using WMS

texts:
I have a shapefile with roughly 6,500 points worldwide in it that I am trying to serve with Geoserver 2.2.1 using WMS. All has been apparently fine with it until I implemented a filtering feature on my client application which uses leaflet. When I append a CQL_FILTER (attribute filter, not spatial) into the WMS query, I noticed missing features when zoomed out. When I zoomed back in, they would sometimes reappear, but not always. See the image below -

At the zoom level on the left, Atlanta is not displayed. When I zoom in, it is. However, sometimes not even the point in Tampa is displayed at the zoom level on the left. If I zoom out 3 more levels, no points at all will be displayed. I'm not sure the issue is the CQL_FILTER param, as with 6,500 points it's hard to notice a few missing points at a global scale, but the particular filter I'm showing here as an example filters down to only 3 features, and when 1 to 3 of them are missing depending on the zoom level it's particularly noticeable.

I can recreate this same behavior in geoserver's own layer preview using the CQL Filter in the advanced map options, so I'm fairly certain it is not my client application.

As for my own steps at debugging this issue, I have tried the default Points style/sld file to rule out my own layer styling. I have turned off all the caching that I am aware of. I double checked that my projections are correct - I created the shapefile in ArcGIS 10 using WGS_1984_Web_Mercator_Auxiliary_Sphere as the projection, and the layer is set to EPSG:3857 in geoserver, which I think is equivalent. I also upgraded from geoserver 2.2 to 2.2.1 and have had the same issue in both. I have also deleted geoserver's spatial index file (.qix) and let it get recreated since I've seen similar issues in Arc with corrupt spatial indexes, but obviously that did not work either.

Here's a snapshot from Geoserver's own layer preview with the CQL Filter on and zoomed into the same area as shown above. The red circle is approximately where I should be seeing another point (Atlanta).

I've tried tweaking every other setting I can think of, but am having no luck. I've also looked up and down the geoserver logs and enabled verbose logging, and do not see any errors/exceptions. I also do not see any bad requests in chrome's developer tools.

If I'm missing any critical info, I'll provide what I can, but this is an internal/non-public application.


texts: Projects feature with attachments from one coordinate system to another

texts:
I'm trying to projects a feature with attachments from one coordinate system to another without loosing the attachments.

I'm using the GeoTagged Photos To Points function to create a point feature class which also contains the photos as attachments.

Unfortunately this function does not provide an option to reproject the output feature class. And if I use the project_management function, I lose the attachments.

Minimal example:

import arcpy
arcpy.env.maintainAttachments = True
# maintainAttachments work for the following function:
# - Append
# - Copy
# - Copy Features
# - Copy Rows
# - Feature Class To Feature Class
# - Merge
# - Table To Table
arcpy.GeoTaggedPhotosToPoints_management(folder, photo_gps, "", "ALL_PHOTOS", "ADD_ATTACHMENTS")
arcpy.Project_management(photo_gps, photo_gps_proj, out_coordinate_system)

But unfortunately photo_gps_proj do not contains the attachments !

How can I project photo_gps and keep the attachments ?


texts: Using 3D Raster Domain to get polygon of raster - but stripping Z and M values

texts:
I'm trying to convert a number of large, floating point rasters into polygons so I can view overlaps etc. I've been using 3D Raster Domain sucessfully to do this. The only problem is that the polygons that are created are very large as a Polygon ZM is created, and I really only need a polygon. Any ideas how to do this? I've tried adding this to my script to disable Z and M values:

#disable M and Z values in output
arcpy.env.outputZflag = "Disabled"
arcpy.env.outputMflag = "Disabled"

... but it didn't reduce the size. I don't want to convert them to integer rasters (the original rasters are too numerous and large in size).


texts: How do I remove the boundaries between municipalities of my shapefile?

texts: I have a shapefile of a province right now. It includes the barangay/village boundaries and municipal boundaries. However, I only need the municipal boundaries. How will I remove the barangay/village per municipality? I tried using the merging feature in the edit menu, but it did not merge those barangays with islands as part. I also tried using the dissolve too and same results appear. What should I do?


texts: QgsRectangle do not completely draw a rectangle

texts:
I create a plugin composed by his own canvas (as widget) and in this canvas the user select a rectangle (for the study zone). I used the RectangleMapTool class.
However, when I draw the rectangle, it is not completed as in the figure.

class RectangleMapTool(QgsMapToolEmitPoint):
    def __init__(self, canvas):
        self.canvas = canvas
        QgsMapToolEmitPoint.__init__(self, self.canvas)
        self.rubberBand = QgsRubberBand(self.canvas, True)
        self.rubberBand.setColor(Qt.red)
        self.rubberBand.setWidth(1)
        self.reset()

    def reset(self):
        self.startPoint = self.endPoint = None
        self.isEmittingPoint = False
        self.rubberBand.reset(True)

    def canvasPressEvent(self, e):
        self.startPoint = self.toMapCoordinates(e.pos())
        self.endPoint = self.startPoint
        self.isEmittingPoint = True
        self.showRect(self.startPoint, self.endPoint)

    def canvasReleaseEvent(self, e):
        self.isEmittingPoint = False
        r = self.rectangle()
        if r is not None:
          print("Rectangle:", r.xMinimum(),
                r.yMinimum(), r.xMaximum(), r.yMaximum()
               )

    def canvasMoveEvent(self, e):
        if not self.isEmittingPoint:
          return

        self.endPoint = self.toMapCoordinates(e.pos())
        self.showRect(self.startPoint, self.endPoint)

    def showRect(self, startPoint, endPoint):
        self.rubberBand.reset()
        if startPoint.x() == endPoint.x() or startPoint.y() == endPoint.y():
          return

        point1 = QgsPointXY(startPoint.x(), startPoint.y())
        point2 = QgsPointXY(startPoint.x(), endPoint.y())
        point3 = QgsPointXY(endPoint.x(), endPoint.y())
        point4 = QgsPointXY(endPoint.x(), startPoint.y())

        self.rubberBand.addPoint(point1, False)
        self.rubberBand.addPoint(point2, False)
        self.rubberBand.addPoint(point3, False)
        self.rubberBand.addPoint(point4, True)    # true to update canvas
        self.rubberBand.show()

    def rectangle(self):
        if self.startPoint is None or self.endPoint is None:
          return None
        elif (self.startPoint.x() == self.endPoint.x() or \
              self.startPoint.y() == self.endPoint.y()):
          return None

          return QgsRectangle(self.startPoint, self.endPoint)

    def deactivate(self):
        QgsMapTool.deactivate(self)
        self.deactivated.emit()

Code where I instantiate the class RectangleToolMap:
canvas = self.dlg.widget
tool = RectangleMapTool(canvas)
tool.rectangle()
canvas.setMapTool(tool)
r = tool.rectangle()


texts: Latitude Longtitude to X, Y coordinate system using python

texts:
I am really new to this stuff and I have to do something for work. I am using python pyproj library to convert lat, long values to x, y coordinates.

myProj = Proj("+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

and tried the following command myProj(50.061665, 19.937332).

Why do I need to specify zone when I create myProj, because the zone should be determinable based on lat long values?


texts: Automating batch load of multiple GPX files into PostGIS?

texts:
I've got 50+ GPX files that I would like to "batch load" into a PostGIS database.  All of the track_points data would be loaded into a "track_points" table (with typical GPS fields like lat, long, elevation, time, etc) and the tracks data would be loaded into a similar, appropriately designed "tracks" line geometry table.

I'd prefer to automate this process so that when I get my next 50+, I won't have to hand enter the data into the database.  I personally prefer to use Python to script these kinds of processes, but any suggestions are welcome.

My general thought process is to:

Get a list of GPX files to process (easy enough through standard Python tools)
Loop through each GPX file and extract/convert the necessary data to PostGIS format
Insert the GPS data into PostGIS using the psycopg Python library

I think I can manage steps 1 & 3, but am wondering if there is a relatively simple method/library that would convert the data (tracks and track_points) into PostGIS format, or simply the tabular form that I could insert into the already created table.

I've already read "Is there a good GPS track analysis library?", "How to build a geographic database of GPS logs?", and "How to extract .gpx data with python", and have looked into the GDAL/OGR and FWTools Python bindings, but don't want to reinvent the wheel is someone already has a good method for this.


texts: Integrating DEM within kriging and other interpolation methods?

texts: I would like to see the difference of including DEM in interpolation from without. Some studies on climate data interpolation suggest that interpolation performs better with DEM. How do I do this in ArcGIS 10.2?


texts: Precision of conversion from WSG84 to ECEF and ECEF to ENU?

texts:
To convert from a global to a local coordinate system, I apply the following transformations.

WSG84 to ECEG:

with

and ECEF to ENU:

Let's assume that the semi-major axes a and e are known with high precision (e.g. a = 6378137.0 and e²=6.69437999014e-3), and the reference position in the second transformation is correct.

What is the precision I can expect from this transformation within a short distance from the reference position (let's say less than a mile)?


texts: Toward a protocol for encoding vector data as image

texts:
This is a follow up to this question: Creating Vector Polygons with rendering performance like GISCloud?

In his answer, Yagi outlines a rationale for encoding geographic information in an image format and decoding it in the browser.  He observes that "Currently, to do this, you have to roll your own".  He also observes that there's currently no standard for this.

Given the awesome performance being demonstrated, it seems like the community might benefit from a standard.  From my understanding of the problem, it sounds like a standard way of dealing with it could be implemented.  Call it B-WFS.

My question, then: what would a useful protocol for encoding vector data as images look like?  Is there something making it too complex to tackle usefully, or is it just a case of "no one's done this yet"?


texts: WMTS from MapProxy in QGIS and/or ArcGIS

texts:
My company wants to use imagery from a another company internally for a contract. That company's WMS has a 5+ second latency. Which I'm told is unacceptably slow. So I set up MapProxy and seeded the cache locally. Latency is now about 1-2 seconds. My boss' boss' boss would like to see if WMTS would be faster. So, I turn on the WMTS layer...and nothing. Researching the problem I see that MapProxy defaults it's origin to the lower left corner, but WMTS expects upper left. So I changed that setting, deleted the cache, and made a new one using the new data. Still no WMTS. QGIS seems to be sending requests and MapProxy is reporting 200, but I see nothing rendered. ArcGIS failed to understand the request at all it seems, at least it wasn't able to read the bounds. I know I'm adding the layer correctly in QGIS (2.0.1) because I can add publicly available WMTS services no problem. When I copy the url into the browser it displays the image no problem. I've tried both Restful and KVP versions. Nothing. MapProxy is running on Centos 6 if it matters.

I have attached the GetCapabilities.xml documents. (Identifying information has to be removed as per company requirements.)

Restful WMTS Service

KVP WMTS Service

MapProxy Configuration:

The WMTS is simply turned on, no special modifications.

Cache:

caches:
  Imagery_Cache:
    meta_size: [5, 5]
    meta_buffer: 20
    format: image/jpeg
    request_format: image/jpeg
    grids: [my_grid]
    sources: [Imagery_source]
    cache:
      type: file
      directory_layout: tc
      directory: /data/mapproxy/folder

sources:

sources:
  Imagery_source:
    type: wms
    supported_srs: ['EPSG:4326']
    req:
      url: "https://is.this/a/url?"
      layers: Imagery
    http:
      ssl_no_cert_checks: True

grids:

grids:
  my_grid:
    base: GLOBAL_GEODETIC
    bbox: [##.# ##.#,##.# ##.#]
    bbox_srs: EPSG:4326
    origin: ul

globals:

globals:
  cache:
    base_dir: './cache_data'
    lock_dir: './cache_data/locks'
    tile_lock_dir: './cache_data/tile_locks'
    meta_size: [4, 4]
    meta_buffer: 80

  image:
    resampling_method: bicubic
    stretch_factor: 1.15
    max_shrink_factor: 4.0
    paletted: false
    formats:
      custom_format:
        format: image/png
        mode: P
        colors: 32
        transparent: true
        encoding_options:
          quantizer: fastoctree
      image/jpeg:
        encoding_options:
          jpeg_quality: 60

So, any thoughts as to what is wrong and how I can get the WMTS to serve up properly?


texts: How to link a pdf file to a specific feature with QGIS?

texts:
I'm trying to link a pdf to a parcel (see image: http://imgur.com/RbauB3v) where you click on one parcel and a specific pdf comes up and if you click on another then a different pdf comes up.

I've got it so far that if you click on a parcel it shows a pdf but the problem is that is shows the same pdf for every parcel!

Is there anyway I can assign it to just one parcel by lets say the objectid?

I'm having trouble setting the action so that it opens each parcel.


texts: How to offset/displace overlapping lines in a web map?

texts:
I have a dataset of points representing a To/From status - something like distribution centres and delivery addresses.

When the user clicks on a point I'd like to show the relationships between this point and its related points using lines. Eg when clicking on an address, show the path to the delivery centre.

This works fine, except when there are multiple overlapping lines. In this situation I'm looking for a method to separate the lines automatically, so that it's clear there are multiple lines, and each can be interrogated.

I've put a simple example on JS Fiddle - this sample creates 5 overlapping lines. Is there any way to displace them automatically, while retaining the start and end nodes?

This is the effect I'm aiming for, ideally handling multiple overlapping lines automatically:

This example uses Google Maps, but an answer in ArcGIS Server's JS API would be acceptable too.


texts: Can I use the LEFT (string) operator in ExecuteSQL()?

texts:
In QGIS I can use the LEFT operator on a dataset like this:

LEFT( "Name2",3) ='031'

What is the syntax for using this expression with ogr in python? Isn't this the right syntax?

a=water.ExecuteSQL("select left('Name2',3)='031' from 'benlux_w_min' ")
'ERROR 1: SQL Expression Parsing Error: syntax error`

A simple where clause works just fine

In [53]: a=water.ExecuteSQL("select * from 'benlux_w_min' where 'Name2'='0310027' ")

In [54]: a.GetFeatureCount()
Out[54]: 367


texts: Seeking OpenJump tool for point/node conflation/matching?

texts:
Are there any OpenJump tools for conflating (matching) two layers of point (node) data?

I'm particularly interested in using such a tool to sync data with OpenStreetMap, where I've created a wiki page on the subject of conflation. That page links to some university research projects on the conflation of vector datasets using OpenJUMP, however they don't share their work.

I know algorithmically this task shouldn't be that difficult especially compared to handling generic vectors, as it would simply use distance and similar properties along with some math, however I suppose I'm more interested in a graphical interface implemented in open source software.

I wrote up a possible workflow for conflating nodes here on the OSM wiki. I thought of trying to implement it in JOSM.

I think it might be useful to say what use cases I'd be interested in using this tool for.

Syncing Virginia interstate exit data from VDOT with OSM
Syncing GNIS feature points with OSM

After an extension to conflate polygons (using their centroids), the tool could also be used for:

Transferring properties (tags, attributes) between nodes representing a house along with addressing data to a polygon (area, closed way) of the actual building

I'm sure there are many other possibilities, but these are the few I'm planning to use  it for.


texts: Can the kml file displayed in Google Earth for Android

texts:
I have created some 3D buildings and route in kml format and it worked properly in Google Earth (desktop). When i put it in Google Cloud Drive and open it by Google Earth for Android, none of the model can be seen.

Grateful if anyone suggested me whether Google Earth for Android support this KML file or Are there any changes in coding is required.

Besides, are there any resources available from Google to test the kml file in Google Earth for Android.


texts: MakeRouteEventLayer_lr Vs Display Route Events

texts:
I've two script tools. The first one prepares data for the second script to work. The first tool routes a table along a line using MakeRouteEventLayer_lr. The second script tool reads the selection made on the layer created by first tool.

First tool:

SMr = arcpy.MakeRouteEventLayer_lr(Route,"ID",dbf,"ID LINE Beg End","Lyr Events")

SMrLyr = SMr.getOutput(0)

SMlyr=arcpy.mapping.AddLayer(df,SMrLyr)

Second tool:

Route  = arcpy.GetParameterAsText(0) # I point to the layer created from First tool

SMSelCount = int(arcpy.GetCount_management(Route).getOutput(0))

If I select few records in the layer from first tool then the second tool should count it for me but instead it gives me count of entire table.

Now if I use Display Route Event(by right clicking the table) and use the resulting layer for second tool then every thing works fine. The second tool can read any  selection made on the layer.

So I was wondering if there is any difference between these two ways of routing. I was under the impression that they are essentially the same thing. Could it be because of the workspace, one using "in-memory" vs other using default.gdb in C drive?


texts:
Using ArcGIS, QGIS, Grass, and/or GVSIG:

What are some of the tools and processes involved in building effective heat maps?
What are the plugins involved?
What are the major data requirements?
What are some of the flaws with existing heat maps?
What are some of the issues that heat maps cannot cover effectively?
How not to do a heat map?
Are there better alternatives (in the same context) than heat map for data representation?


texts: Using gdal_translate: how to extract a raster layer into a new tif and keep the original's metadata

texts:
I am trying to extract each raster layer from a geoTiff into a separate tiff. I can do this quite easily by doing

gdal_translate -b 1 input.tif output.tif

However I need to preserve the metadata from the original geoTiff into the separated raster tiff. I have been playing around with the API on this link (http://www.gdal.org/gdal_translate.html) however have been unsuccessful

Any idea on how to do this?


texts: How to automatically format the legend labels? in QGIS

texts:
I have to do 10 maps (for a start) with "style" -> "graduate" with quantile mode.

I need to format the label of the different ranges so that the legend looks comprehensible:
Instead of having 1.0000 - 7.0000 I would like to have 1% - 7% without me typing.

In arcgis there is the possibility to format the labels. But I cannot seem to find anything similar in QGIS 1.8. Is it possible with some plugin or other QGIS version ?

Update some years later: the new version seems to have fixed the problem as the answer belows points out.


texts: Using python and ArcPy to read multiple files, apply a function and write into new a directry

texts:
I haven't used python yet; I want to read multiple GeoTiff files and resample them based on another single file, and be able to write all the files in a new directory. I can use R raster package to that, but it takes a long time. Here is what I do in R and I want to be able to do the exact process in Python using ArcPy module.

    setwd("path to directory")
    files <- list.files(path=".", pattern="tif", all.files=FALSE, full.names=TRUE)
    s <- stack(files)
    resamp <- raster("new_raster.tif")
    s2 <- resample(s, resamp,method='bilinear')
    s3 <- unstack(s2)
    setwd("new_directory)
    for(i in seq_along(s3)){writeRaster(s3[[i]], file=outputnames[i],format="GTiff", options=c("INTERLEAVE=BAND","COMPRESS=NONE"))}


texts:
My lat/long points have an offset (2-3 meters) in both X and Y - about 2m SE in my region (-120.644,35.260) - from where the points used to be (unless my memory is failing.) I'm using Google Earth 7.1.

They (exactly lat/long points in GE) plot correctly on an OSM and other GE overlays. OSM shows the point in the correct place on the image but in GE the point appears offset to the SE. OSM streets also have the same offset.

All the historical layers are offset as well.

Has something happened in the last few months to caused GE to plot points incorrectly?


texts: Abort pending WMS image requests on map pan

texts:
Everytime there's a move in the map the OL3 requests a new image to the GeoServer ignoring the previous requests it has done which may lead to multiple pending image requests that are not going to be used.

Is there a way to abort these requests like in jQuery with .abort() for AJAX requests?