texts:
I am trying to write a script that accomplishes 2 tasks: change the extent of all map frames using a boundary layer as reference and rotate the map frames by -17 degrees. To complete the first task, I decided to use the getLayerExtent() method to retrieve the extent of the boundary layer, then assign each mapframes' extent by using setExtent(). Unfortunately, I cannot get past this first task.

When trying to run the getLayerExtent() method in my script I am met with an error that reads:

  ValueError: MappingLayerObject object at 0x000002210CE738A0

Needless to say, the error message isn't very descriptive and doesn't offer many hints on how to tackle it. When searching online, I could not find many examples on how to use this method or even how to fix the value error.

Below is my script where I try to use the specified method:

import arcpy

aprx = arcpy.mp.ArcGISProject("C:/JHU406Lab7/Lab7/Toronto/Toronto_Overview.aprx")

region_map = aprx.listMaps("Region")[0]

boundary = region_map.listLayers("boundary")[0]

lyt = aprx.listLayouts()[0]

mfs = lyt.listElements("MAPFRAME_ELEMENT")

for mf in mfs:

    mf.camera.setExtent(mf.getLayerExtent(boundary))

    print("Map frame extents have been updated")

    mf.elementRotation = -17

    print("Map frames have been rotated")

aprx.saveACopy("C:/JHU406Lab7/Lab7/Toronto/Toronto_Overview_q5.aprx")

del aprx

To use the getLayerElements() method, I first had to retrieve the boundary layer by specifying the map and layer. I then retrieved the layout and the list of map elements I want to iterate. Lastly, I tried to  iterate through each map frame and use the getLayerExtent() method. I'm stumped on what is causing this error.


texts: extract roads, railways and canals out of old orthophoto

texts:
how is it possible to extract the roads, canals and or railways out of this image (black-white)?

Or how is possible to create lines at the edges of the roads, canals and or railways?


texts: Compute area in square meters rather than degrees using PyQGIS?

texts:
If you use the field calculator in Qgis considers the area in square meters, and when using the script gives the area in degrees.
My layer is in the coordinate system WGS 84

vlayer=r"K:/temp/gpp_utm.shp"
vlayer.startEditing()
myField=QgsField('myattr',QVariant.Double)
vlayer.dataProvider().addAttributes([myField])
vlayer.updateFields()
idx=vlayer.fieldNameIndex('myattr')
# vlayer.addAttribute(myField)
# idx=vlayer.fieldNameIndex('myattr')
e=QgsExpression('$area')
e.prepare(vlayer.pendingFields())
for f in vlayer.getFeatures():
    f[idx]=e.evaluate(f)
    vlayer.updateFeature(f)
vlayer.commitChanges()

How to calculate the area in meters and not in degree in my script?
I tried to put $area(m2) but it did not help.


texts: PIP no longer works with QGIS 2.18.x Python (TLS version no longer supported)

texts:
We've run into a problem installing Python packages (using pip) for use with QGIS 2.18.x (error messages copied below)

Our diagnosis is that:

Pypi now requires TLSv1.2 support (http using TLSv1.0 and TLSv1.1 are now blocked). See https://pyfound.blogspot.com.au/2017/01/time-to-upgrade-your-python-tls-v12.html
QGIS 2.18.x Python is 2.7.5 and does not support TLSv1.2, and seems it cannot be easily upgraded

Is there a way to enable TLSv1.2 support in the QGIS 2.18.x Python environment, to enable PIP install to work?

PIP install error message from OSGEO4W Shell window
(using pint as example: all packages give similar results:)

C:\Windows\System32>python -m pip install pint
Collecting pint
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\pip\_vendor\requests\pack
ages\urllib3\util\ssl_.py:318: SNIMissingWarning: An HTTPS request has been made, but the SNI (Subject Name Indication) extension to TLS is not available on this platform. This may cause the server to present an incorrect TLS certificate, which can cause validation failures. You can upgrade to a newer version of Python to solve this. For more information, see https://urllib3.readthedocs.io/en/latest/security.html#snimissingwarning.
SNIMissingWarning
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\pip\_vendor\requests\packages\urllib3\util\ssl_.py:122: InsecurePlatformWarning: A true SSLContext object is not available. This prevents urllib3 from configuring SSL appropriately and may cause certain SSL connections to fail. You can upgrade to a newer version of  Python to solve this. For more information, see https://urllib3.readthedocs.io/en/latest/security.html#insecureplatformwarning.
InsecurePlatformWarning
Could not fetch URL https://pypi.python.org/simple/pint/: There was a problem confirming the ssl certificate: [Errno 1] _ssl.c:504: error:1407742E:SSL routines:SSL23_GET_SERVER_HELLO:tlsv1 alert protocol version - skipping  Could not find a version that satisfies the requirement pint (from versions: )

No matching distribution found for pint


texts: How to derive GetFeature URL from capabilities document?

texts:
Given this capabilities definition how can I compose a GetFeature request that downloads the actual geodata? I came up with the following URL referring to WFS 2.0.0:

http://fbinter.stadt-berlin.de/fb/wfs/geometry/senstadt/re_umweltzone2007? \
  request=GetFeature&service=WFS&version=2.0.0&typeNames=fis:re_umweltzone2007&outputFormat=GML2&

For some reason it does not accept &srsName=EPSG:4326 as a parameter. Is it possible to create a download URL such as &outputFormat=shape-zip which I saw elsewhere? By now I use wget to trigger the download.


texts: PostGIS - identify duplicate points from polygons/linestrings?

texts:
I'm trying to identify polygons/linestrings with duplicate points, using PostGIS.
I just want to tell whether the geometry has duplicate points, not remove the extra point.
I know there is ST_RemoveRepeatedpPoints, but as far as I can tell it only does something if the points are right after one another i.e.
ST_RemoveRepeatedPoints(st_geomfromtext('Polygon((0 0, 0 1,
                    0.5 1, --> 3. point
                    0.45 0.45, 0.5 0.5,
                    0.5000000001 1, --> notice it is basically the same as 3. point
                    1 1,  1 0, 0 0))'));
-- Answer:   'Polygon((0 0, 0 1, 0.5 1, 0.45 0.45, 0.5 0.5,  0.5000000001 1,1 1,  1 0, 0 0))'
-- nothing changed

-- same before, but with added repeating point
ST_RemoveRepeatedPoints(st_GeomFromText( 'Polygon((0 0, 0 1, 0.5 1, 0.45 0.45, 0.5 0.5,  0.5000000001 1,1 1, 1 1,  1 0, 0 0))'));

-- Answer:   'Polygon((0 0, 0 1, 0.5 1, 0.45 0.45, 0.5 0.5,  0.5000000001 1,1 1,  1 0, 0 0))'
-- repeating point (1,1) is removed, but nothing else changed

otherwise the function could be used together with ST_NPoints to see if there's a difference between the original and modified geometry like this:
ST_NPoints(ST_RemoveRepeatedPoints(geom)) <> ST_NPoints(geom)

I have also tried ST_SnapToGrid to identify geometries with duplicates, but it seems to have the same behavior as ST_RemoveRepeatedPoints.
So far the simplest solution I have found is to

dump the geometry to points
collect the points and then
use ST_SnapToGrid to see if any points are removed, i.e.:

WITH dumppts AS (SELECT (ST_DumpPoints(geom)).geom FROM some_table)
, collect_pts AS (SELECT ST_Collect(geom) AS geom FROM dumppts)

SELECT
ST_NPoints(st_removerepeatedpoints(geom)) <> ST_NPoints(geom) -- returns TRUE if there are duplicate points
FROM collect_pts

but this seems rather slow, because you have to dump and then collect the geometry back together.
Is there a quicker/smarter solution?
EDITED: modified the example geometry so its now polygons


texts: Shapefile and geocoding to retrieve an area name

texts:
I'm very new to GIS and all things geospatial!

I have dbf, prj, sbn, sbx, shp and shx files that have been provided to me. The shapefile contains information on parts of the country that have been segmented into logical areas. These are borders which do not overlap.

I want to take a geocoded location and some how determine which area the location falls within.

I'm unsure as to how to proceed; are there libraries which directly read shapefiles, should I convert to KML, or to a database.

Platform is PHP.

If you require any more information please let me know.

EDIT Sorry I had no idea how much info you guys needed so I'll update here based on comments. No, I have no software at all, I was just provided these files and asked to "make it work". Never done any geospatial work before.

EDIT2 I downloaded QGIS and managed to open the shapefile and I can see my state. I'll explain precisely what I am trying to do: the shapefile has the state split into various administrative areas. What I need to do is to take the geocoded address that someone enters and tell them what area they live in.

Here is the datastructure they have provided me:

Formats available: MAPINFO TAB; ESRI SHAPE
Datum: GDA94
Projection: MGA94 Zone 55
Geometry: Polygon
Field Type Description
ID INT32
NAME CHAR(60)
PLAN_REF CHAR(8)
COMMENTS CHAR(50)
UFI CHAR(12)
FMP CHAR(12)
CREATED_ON CHAR(19)


texts: Extracting values From several raster .tif files to a .csv file using python

texts:
I have multiple raster .tif files. I am using Python, and I want to extract cell's attribute values from all rasters based on the latitude and longitude similarity, and saved in a .csv file, the output should be formatted as follows:

lat lon value* value**

20  77  12     20

*denotes attribute value from raster 1

**denotes attribute value from raster 2

Can you please help me in this ?


texts: Export image collection to pixel value across bands

texts:
I am trying to export an image collection to a 2D table where :

axis 0 are pixels over the collection
axis 1 are bands from the image collection

I would loose spatial information, as pixels would be unordered which is OK. It could look like this (pixel Id unecessary):

I found the function toArray() that seems to do exactly what I want, but I am unable to export it to a CSV file or any type of format

// Load a Landsat collection, select the bands of interest.
var imageCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .select(['B1', 'B2', 'B3'])
  .filterDate('2017-09-01', '2017-09-02');

// Make an Array Image, with a 1-D Array per pixel.
var im_array = imageCol.toArray();

I have tried to transform this into a featureCollection but it doesn't work :

var im_feats = im_array.reduceToVectors()

I would have then export it in CSV with Export.table (it would be quite a heavy file):

Export.table.toDrive({
  collection: im_feats ,
  description:'pix_bands',
  fileFormat: 'CSV'
});

Am I taking a wrong approach ? Is there a trick to get around this issue ?


texts: Using mixed fonts for one text field when publishing to ArcGIS for Server?

texts:
I have the following name field, which appears perfect in ArcMap when i use an Armenian font.

Նոր Նորքի 9-րդ զանգված

The problem is when i publish this as a service the Latin characters appear as Նոր Նորքի □□րդ զանգված when using Armenian or the when using Latin font the armenian text appears □□ □□□□□ 91□□ □□□□□□□

Is there any way of being able to fix this?  I believe i need to use 2 different fonts


texts:
I am a pretty new user to Earth Engine, and I am trying to map over a list in order to use the numbers inside as indices. As a quick example, in R this would be

for (i in 1:10){print i}

In Earth Engine, doing this separately works beautifully. Here, "points" is a FeatureCollection made up of 10 points, and my goal is to filter the "l8sr" ImageCollection to only give me the images that pertain to point #2. Then, I get mean NDVI for the time series.

var coords = points.geometry().coordinates()
  print('coords', coords)

  var onepoint = ee.Geometry.Point(coords.get(2)) //for example
  print(onepoint)

  var test = l8sr.filterBounds(onepoint)
  print('test', test)
  print(ui.Chart.image.series(test, onepoint, ee.Reducer.mean(), 30))

BUT, when I try to map this over all 10 points in "coords", I get a list of 10 empty ImageCollections (I can't print to a chart because these are undefined). Why is this?

Here is reproducible code:

// Make a list of Features.
var features = [
  ee.Feature(ee.Geometry.Point(97.58, 24.36), {name: 'sample1'}),
  ee.Feature(ee.Geometry.Point(99.45, 20.26), {name: 'sample2'}),
  ee.Feature(ee.Geometry.Point(93.24, 21.42), {name: 'sample3'}),
  ee.Feature(ee.Geometry.Point(96.13, 24.12), {name: 'sample4'}),
  ee.Feature(ee.Geometry.Point(96.88, 22.48), {name: 'sample5'})
];

// Create a FeatureCollection from the list and print it.
var points = ee.FeatureCollection(features);
print(points);

var l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2015-06-23', '2020-01-31')
                  .filterBounds(points)
print("Landsat", l8sr)

//this works (down to /$/$)
var coords = points.geometry().coordinates()
print('coords', coords)

var onepoint = ee.Geometry.Point(coords.get(2))
print(onepoint)

var newt = l8sr.filterBounds(onepoint)
print('newt', newt)
print(ui.Chart.image.series(newt, onepoint, ee.Reducer.mean(), 30))

// /$/$ /////////////

// this doesn't work
var ser = ee.List.sequence(0,4)
print(ser)

var test1 = ser.map(function(q){
  var sub = l8sr.filterBounds(ee.Geometry.Point(coords.get(q)))
  return sub
  //return ui.Chart.image.series(sub, point_a, ee.Reducer.mean(), 30)
})

// test1 is a list of ImageCollections, but each one has no actual bands.
print(test1)


texts: Reprojecting basemap layers in ArcMap without slowing down zoom and pan options?

texts:
I am using ArcGIS 10.5.1 on a high processing computer. I have 4768 data points that I am required to zoom into and pan around to determine their composition.

They are projected in NAD 83 CSRS Statistics Canada Lambert and must stay in this projection.

When I reproject the basemap to this projection, it significantly slows down the ability to zoom in and out and pan around.

Is it this projection that is the issue, or the order in which I am doing things?


texts: Selecting features from one layer into new layer based on IDs using PyQGIS?

texts:
I'm completely new to PYQGIS. I'm pretty handy with MapBasic in MapInfo and was wondering how I would 'translate' this command from SQL select into Python (hope that makes sense)

Essentially, I have a layer with multiple objects - but I only want to select and extract the objects which have an id like %ABC% into a separate layer.

In MapBasic I would:

Select * from *examplelayer* where id Like "%ABC%" into *ABC*

where 'id' is a field I want to search in the table/layer 'examplelayer', id Like "%ABC%" meaning that any id's containing ABC will be selected into a new table/map layer named "ABC"


texts:
Techniques to map small water bodies in desert

I have task to map water bodies in a desert. I used landsat 8 and calculated three indices called NDMI, Tasseled Cap Wetness, & AWEIsh. I am unable to map outlines water bodies as those appear in satellite view of GEE. Screenshot, kml file and code link attached. I used threshold values from literature even Ostu segmentation to find optimal threshold.

What I should do to get map water bodies along with their natural boundaries? Unmixing? High spatial resolution data? Other image classification technique? use https://imagej.nih.gov/ij/plugins/versatile-wand-tool/index.html
https://code.earthengine.google.com/93ceae78cac565d6f8975da1aff1c370


texts: When validating geometries, what are "Nested shells"?

texts:
I am validating shapefiles polygon geometry using shapely, and have also attempted it using FME. Some polygons return the error "Nested Shells".

I cannot find a definition of what the problem is, whether in FME documentation (here) or in google search results on shapely or OGC geometry validation.

Would anyone have a definition or (ideally!) a graphic example?


texts: How to emulate a click on a MapBox GL JS map?

texts:
I need to emulate a user click on my MapBox GL JS map.

I doesn't matter where in my map: could be useful to do in the point at the center of my maps if this could be easier.

UPDATE

I've found here https://stackoverflow.com/questions/18323105/programmatically-opening-a-popup-in-a-mapbox-map something that could be useful ...

map.fireEvent('click', {latlng: L.latLng(28.04419, -81.947864)});

... but seems to be something about Leaflet

I'd like my code coud be something like this ("broken code" ... )

    bounds = map.getBounds();
    x = bounds._ne.lng;
    y = bounds._ne.lat;

    map.fireEvent('click', {latlng: L.latLng(x, y)});

Does fireEvent exist in MapBox GL JS or are there alternatives?


texts: How do you deal with (repair/avoid) broken MXD\LYR paths?

texts:
In relation to this question ...

In the past,  whenever I cleaned up a bit and moved layers around (or changed their names) my MXD's paths were broken and it was time-consuming to correct them.

I'd like to hear from your experiences and maybe learn a few tips from you:

How do you deal with broken MXD paths?

(I'm looking for solutions for one computer or a simple network drive, so if your solution deals with larger networks or the web please say so).


texts: Create PostgreSQL table from GeoJSON data FeatureCollection

texts:
I would like to store GeoJSON data into PostgreSQL database.

WITH data AS (SELECT '{ "type": "FeatureCollection",
    "features": [
      { "type": "Feature",
        "geometry": {"type": "Point", "coordinates": [102.0, 0.5]},
        "properties": {"prop0": "value0"}
        },
      { "type": "Feature",
        "geometry": {
          "type": "LineString",
          "coordinates": [
            [102.0, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]
            ]
          },
        "properties": {
          "prop0": "value0",
          "prop1": 0.0
          }
        },
      { "type": "Feature",
         "geometry": {
           "type": "Polygon",
           "coordinates": [
             [ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0],
               [100.0, 1.0], [100.0, 0.0] ]
             ]
         },
         "properties": {
           "prop0": "value0",
           "prop1": {"this": "that"}
           }
         }
       ]
     }'::json AS fc)

SELECT
  row_number() OVER () AS gid,
  ST_AsText(ST_GeomFromGeoJSON(feat->>'geometry')) AS geom,
  feat->'properties' AS properties
FROM (
  SELECT json_array_elements(fc->'features') AS feat
  FROM data
) AS f;

I can query but I need to save data as a table.
Does anyone know how to store this kind data into a table in PostgreSQL.?

Reference to
Storing GeoJSON FeatureCollection to PostgreSQL with PostGIS?

This reference does only query but I need save the layer as a table.


texts: Convert a WFS GET request to a WFS POST request

texts: I have a WFS GET Request, created by A WFS Layer in QGIS. However, for my program I need to have a POST Request. How do I convert the WFS GET request to a WFS POST request?


texts: Compiling a list of units that border each other

texts:
I have a shapefile of U.S. counties.  I am interested in compiling a list of those that border each other.

In particular, I am interested in those that border each other but are in different states.  For example, I am interested in knowing that Rabun County borders counties in North Carolina and which ones and not the counties in Georgia that it may also border.

I work with ArcGIS 10.1 and use ArcMap.

I see that there is a feature called Polygon Neighbors.  However, I don't fully understand how to use it.  I'm not sure what to put into the "Input Features" or how to proceed from here.

I would very much appreciate any help!


texts: Conversion from number field to text field

texts:
I have a field "folio" formatted as a number with decimal. I added a new string field "folio_1" and using field calculator to populate it with the value of "folio". The problem is that the new field drops the zero at the end of all the numbers. How can I conserve the integrity of the number with decimal when converting a number field to a string one?

Folio            Folio_1
32901.0640       32901.064


texts: Arcade expression to display alias field name and value in popup

texts:
I'm trying to return certain fields in my popup in AGOL depending on the value in the visit type field. I'm struggling to figure out how to return the alias field name and value. Right now my expression only returns the field name. Any suggestions?
var list1 = ['COMNAME','ABUNDANCE']
var list2 = ['CLASS']
for (var i in $feature) {
    if ($feature.VISITTYPE == 'Mapping') {
        return list1;
    }
    else if($feature.VISITTYPE == 'Treatment'){
        return list2
    }
}

Current results in popup:
 - COMNAME
 - ABUNDANCE

If I change list1 to :
var list1 = [DomainName($feature,"COMNAME"), DomainName($feature,"ABUNDANCE")]

Then I get the values now but without the field alias name:
- Black mustard
- Linearly Scattered

My desired results in popup is below:
   Common Name: Black mustard
   Abundance: Linearly Scattered


texts: Where does Open-Streetmap forest data come from?

texts:
If you open this url

http://www.openstreetmap.org/?mlat=50.6524658203125&mlon=9.4482421875&zoom=13

you are in a forest in germany. I downloaded the complete Open-Streetmap Data for Germany. I can find the waynodes for the tracks in this area... But not for the forest itself

Cna anyone help? Which way (id) defines the forest under the marker?

Or is it just the Background image for areas where's nothing defined?


texts: Why does QGIS Polygonize ignore the no data value?

texts:
I'm running into an issue georeferencing a number of PNG image files. The files simply contain an image where each pixel is either transparent or a single color, forming a coverage mask. In the past, when I've georeferenced these images, the 'use 0 for transparency' checkbox masked off the transparent parts of the images and I was able to polygonize the colored parts. However, with a number of images in a recent patch QGIS will simply not mask off the transparent sections. I've opened the files with an image editor and resaved them to no effect, so it's not a formatting problem with the input png.

I can somewhat work around this by setting the value of the transparent section to the no data value, thus masking off the image in the raster layer. However, when I polygonize the raster, the no data value still gets polygonized into the .shp file.

How do I polygonize the raster and specify only the single band I want to polygonize?


texts: Efficiently load pgpointcloud from PostGIS table values

texts:
Is it possible to efficiently load data into a pgpointcloud PCPatch table from an xyz type of table within the same database?

The github read me says - “Usually you will only be creating tables for storing PcPatch objects, and using PcPoint objects as transitional objects for filtering...”

And

“PostgreSQL Pointcloud storage of LIDAR works best when each "patch" of points consists of points that are close together, and when most patches do not overlap.”

So does this mean I should convert each individual xyz record to a PCPoint record, then convert n records to a patch? Window function?

Or should I create an array from the xyz records, then create a patch from the array?

I’m not sure what the most efficient/effective way to undertake this and ensure that the patches are made up of close together points.


texts: problem to Disable the security in server.dat file

texts:
I am trying to remove the security from web services.While seaching , I found that we have to set

<SecurityEnabled>True</SecurityEnabled>

                  To

<SecurityEnabled>False</SecurityEnabled>

Now after replacing this true with false , I saved the file. Now when I restart my computer and reopen this file, It is again set to <SecurityEnabled>True</SecurityEnabled> .I am unable to understand why this is happening?


texts:
I use WFS=1.0.0 on the Mapserver v5.6.6 and v6.0.1 with MS4W package. I want to get the Intersects of Spatial Operator. When I make getcapabilities request I see just "Intersect" like:
ogc:Filter_Capabilities>
-
-

When I change version of WFS to 1.1.0 I see the Intersects Spatial Operator.
I need get the Intersects on WFS=1.0.0. Is it possible? If so how can I get it?
P.S. When I use Geoserver then it seems to work.


texts: Creating points on either side of the line at specified intervals

texts:
With the tool "Generate points along line", it creates points along the line which I don't want.

I Want the distance between two alternating points (zig zag) should be 250m and between the distance from point to line should be 50m.


texts: Trying to completely Exit from the MapInfo application?

texts:
I tried the following two things:

 include "mapbasic.def"

 Declare Sub Main
 Declare Sub Exit

Sub Main()

    If WindowInfo(WIN_MAPBASIC, WIN_INFO_OPEN) then
        Close Window WIN_MAPBASIC
    End If

 End Sub

and

 Sub Exit()
    End Program
 End Sub

But my MapBasic program does not completely exit from MapInfo (close MapInfo).


texts:
I have a NDVI raster file generate from Landsat-8 data and using the Identify Features tool I am able to observe the individual values per pixel. I also have a vector kml file  breaking down a region into a number of areas. I am trying to use zonal statistics to gain average, mean etc from the raster data for each area of the vector file. However, when I complete the zonal statistics all the calculated values are null.
Any ideas as to what I might be doing wrong?


texts: Replace a feature class in a feature class list

texts:
I'm using arcpy with arcgis 10.0. I try a script that erase successively all the feature class of a list depending on a reliability field value. Each time a feature class is erased by an other one, I want to replace the corresponding feature by the erased one. The problem is that I'm doing it wrong. Check the "HERE IS THE PROBLEM" comment in the code below :

import arcpy, os
arcpy.env.workspace = r'Z:\Documents\SIG\tests\agregation_auto\erase.gdb'
fcList = arcpy.ListFeatureClasses()
outFolder = r"Z:\Documents\SIG\tests\agregation_auto"
count = 0 #help select the a feature class
cnt = 1 #help select the b feature class
tList = len(fcList) #Count the number of elements in the list
for count in range(tList - 1):
    for cnt in range(cnt,tList):
        desca = arcpy.Describe(fcList[count])
        descb = arcpy.Describe(fcList[cnt])
        namea = desca.baseName # Get the name of the a feature to name the out_poly
        nameb = descb.baseName # Get the name of the b feature to name the out_poly
        arcpy.MakeFeatureLayer_management(os.path.join(arcpy.env.workspace, fcList[count]), "fca")
        arcpy.MakeFeatureLayer_management(os.path.join(arcpy.env.workspace, fcList[cnt]), "fcb")
        arcpy.SelectLayerByAttribute_management("fca","NEW_SELECTION")
        arcpy.SelectLayerByLocation_management("fcb","INTERSECT","fca")

        if arcpy.Describe("fcb").FIDSet: #Check for selection
            rowa = arcpy.SearchCursor("fca")
            rowb = arcpy.SearchCursor("fcb")
            for row in rowa:
                fiabia = row.fiabilité #get the int value of reliability (all entities of a layer have the same)
            for rows in rowb:
                fiabib = rows.fiabilité
         if fiabia > fiabib:
             out_poly = os.path.join(outFolder, nameb + "_er_" + namea + str(count) + "_" + str(cnt))
             tmpb = arcpy.Erase_analysis("fcb","fca", out_poly)
             fcList[cnt] = tmpb # I think HERE IS THE PROBLEM
          if fiabia < fiabib:
             out_poly = os.path.join(outFolder, namea + "_er_" + nameb + str(count) + "_" + str(cnt))
             tmpa = arcpy.Erase_analysis("fcb","fca", out_poly)
             fcList[count] = tmpa # I think HERE IS THE PROBLEM


texts: Make shapefile in ArcMap Editable in ArcGIS Collector

texts:
I created a shapefile called random points in ArcMap. I assigned it a coordinate system, and I then shared as a hosted map service. Under Feature Access I checked Create, Delete, Query, Sync, and Update. Under the Sharing tab I checked the appropriate group.

When I accessed the Collector App, I was able to see the points I added as well as see its proper spatial location.

Though I was able to delete or edit a point, I was not able to see the add point feature. Under settings I also checked Add, Update, and delete features.

When I add another feature class I see the Add point feature in the collector app, but for some reason my recent created layers don't have the tab shown up.

I wonder if it's something I did not check off or took into consideration when creating my layer.


texts: Adding points defined by Latitude and Longitude to map using QGIS?

texts:
I have an OS map on which I want to display points defined by Latitude and Longitude, but the map coo-ordinates have values I do not recognise or understand. The image shows the map, with my first point with its latitude and longitude. It is shown where I would expect it to be, but it does not appear on the map.. At the bottom of the image the co-ordinate system of the map is something I do not recognise, certainly not latitude and longitude. Please, how do I set it up so that my points are correctly displayed?

I am using a csv file of latitudes and longitudes, the first few lines are

Latitude,Longitude
50.82932,-1.05779
50.82934,-1.05790
50.82933,-1.05791
...

The base map is from the UK Ordnance Survey, starting at: this entry point, downloading the data for the National Grid Reference square SU, and selecting the map for su60se.

I tried to add the file of data to the map using "Add Delimited Text Layer." I did check, using Google Maps, that the first co-ordinates do represent the point shown in my image. I have set the map Co-ordinate Reference System, CRS, to WGS84. The co-ordinates in the red ellipse appear to be with reference to the top left corner of the map as 0,0 (possibly pixels?)

I should also mention that this is my first attempt to use QGIS, so please forgive my fumbling approach.


texts: Creating multiple users for postgis db created with qgiscloud plugin?

texts:
I have created a PostGIS database using qgiscloud plugin.

I would like to give access to 3 more users (usernames and associated passwords) in order for them to update 3 specific layers only, without having privileges to edit other layers.

Is it possible and how?


texts: Labelling line layer using QgsPalLayerSettings() QGIS failed

texts:
I'm trying to set label on polygon and line layer using PyQGIS. Setting label on polygon layer succeed, but fails on line layer.

Here my script :

        #Set polygon layer
        palyr = QgsPalLayerSettings()
        palyr.readFromLayer(layer_base)
        palyr.enabled = True  # this works
        palyr.fieldName = 'no_urut'  # this works
        palyr.fontSizeInMapUnits = False
        palyr.textFont.setPointSize(6)  # results in 4 - seems to be integer only
        palyr.textColor = QColor(0, 0, 0)  # this works
        palyr.writeToLayer(layer_base)

        #set label line layer
        palyr3 = QgsPalLayerSettings()
        palyr3.readFromLayer(layer_jalan)
        palyr3.enabled = True  # this works
        palyr3.fieldName = 'nama_jalan'  # this works
        palyr3.fontSizeInMapUnits = False
        palyr3.textFont.setPointSize(11)  # results in 4 - seems to be integer only
        palyr3.textColor = QColor(0, 0, 255)  # this works
        palyr3.placementFlags = QgsPalLayerSettings.AboveLine
        palyr3.writeToLayer(layer_jalan)

Any ideas ?, btw when I load it and check on the layer properties, the label has been set according to the properties in the script above and, when I push the apply button on this properties window, the label appears on the canvas.


texts: QGIS - processing area and calculating volume of polygon with know elevation from a DEM

texts:
I got a DEM of a watershed and I need to create a dam at a fixed elevation (with a rectangular polygon corresponding at the top of the dam) and calculate it's area (including the slope (fixed value) between the perimeter of the top of the dam and the DEM) and the corresponding volume.

Do you have any idea of the tools I can use for :

processing the real area of the dam (ie : including the area of the slope)?
calculating the volume of the full dam


texts: python scripts for accessing info from a postgresql database in case of a qgis application?

texts: I am beginner and I need to know if I could access a postgresql database from python script made for a qgis application


texts: How to get Leaflet interactivity from Geoserver WMS?

texts:
I have geoserver serving shapefiles as a WMS into Leaflet.

I'd like to add interativity that shows a popup with WMS attributes in the Leaflet map when the user hovers over or clicks on a WMS feature.

I've seen the wax/leaflet interactivity, but is that possible with WMS and Leaflet?


texts: Updating a table using outputs and insert.da.cursor

texts:
I have the following code however I need the outputs in a table which can be analysed further. How can I take the outputs from the code below and input into the table.
Ive used insert cursor but am getting the following error:

Runtime error
Traceback (most recent call last):
  File "< string >", line 11, in
TypeError: value #1 - unsupported type: Result

import arcpy
arcpy.env.workspace = r"D:\Zahir\AIS\Raw_AIS_Data\2013\Test\\"
gdblist = arcpy.ListWorkspaces()
table = arcpy.CreateTable_management("D:\Zahir\AIS\Raw_AIS_Data\Test\\", "Table.dbf")
arcpy.AddField_management(table, "Name", "TEXT")
arcpy.AddField_management(table, "PosCount", "SHORT")
arcpy.AddField_management(table, "TransCount", "SHORT")
fieldlist = ["Name", "PosCount", "TransCount"]

for g in gdblist:
    desc = arcpy.Describe(g)
    name  = desc.basename
    print name
    trans_table = "main.tblTransit"
    pos_table = "main.tblPositions"
    position_count = arcpy.GetCount_management(g + "\\" + pos_table)
    transits = arcpy.Frequency_analysis(g + "\\" + trans_table, "Trans_freq", "MMSI_ID")
    transit_count = arcpy.GetCount_management(transits)
    print g + " Unique Positions = " + str(position_count) + " and Vessel Count = " + str(transit_count)
    with arcpy.da.InsertCursor(table, fieldlist) as cursor:
        cursor.insertRow((name, position_count, transit_count))


texts:
I use QGIS 2.18.5 and I get in stuck with two polygons I cannot solve.
The two polygons I have looks like:

What I want to get is circled with red:

What I tried/made: Cut 1 polygon and paste to other layer. Used polygon to lines toolbox then explode lines, then deleted the lines which do not needed then used the add feature to connect lines. Last I dissolve new lines then convert lines to polygons. The new polygon looks like that I wanted to get in QGIS but when I use wayline converter which convert polygons to my tractor GPS it is only shows that polygon what was originally on that layer when I started.
Look I deleted the line I do not need:

When I select Vertex Tool it appears the line I deleted earlier:

I think this is the main problem but I do not know how to remote it well..
How can I solve this problem?


texts: CartoDB: Dynamic text with different links in infowindows

texts:
I have a set of text with links for a set of countries that I'd like to display on click in CartoDB. I've created a column "text" and it displays fine on click. The problem is that I can't show the links at all.

Is there a way to keep them active? The a href tag isn't working since the links are always at a different position within the text.


texts: How to send data parameter on leaflet-ajax

texts:
i have a problem about sending data parameter to leaflet ajax,
i'm using https://github.com/calvinmetcalf/leaflet-ajax
anyone know how to send data parameter ? i have reading documentation but did't get it if i using https://github.com/calvinmetcalf/leaflet-ajax

here's my code:

function popUp(f,l){
var out = [];
if (f.properties){
    for(key in f.properties){
        out.push(key+": "+f.properties[key]);
    }
    l.bindPopup(out.join("<br />"));
}}
var myurl = "127.0.0.1/data.json";
var jsonTest = new L.GeoJSON.AJAX(myurl,{onEachFeature:popUp}).addTo(m);


texts: How to check my EPSG and units in R when it is not defined in proj4string?

texts:
After reading my .shp in R using readOGR(), I tried to figure out what is my EPSG and what are my units in my projections?

when I tape

my.shp = readOGR(dsn=setwd(dsn = dsn,
                     my.shp)

proj4string(my.shp)

the returning information is:

[+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0]

so NO information about EPSG or units? are they in miles or in kilometers?

Moreover, when checkinhg my data with summary(my.shp)
my data seems to not be projected (Is projected: FALSE), but they still have proj4string?

Is projected: FALSE
proj4string :
  [+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0]

Thus: are my data projected or not? and how can I find out what is my EPSG and my map units? I need to have my map in Miles, so do I need to transform spTransform() them? If yes, how?


texts: Parallel process for spatial weights matrix production in R

texts:
I have an sf object pData that is rather large (300,000 polygons) and I am trying to produce spatial weights for each polygon:
tilesNb=knn2nb(knearneigh(st_centroid(pData), k = 2))
tilesWeights=nb2listw(tilesNb, style="W")

But the first line ran for hours and then timed out. Is there a way to conduct this step through a parallel process? I tried modifying code from here, but I couldn't quite figure out how to modify their technique for my simpler approach to neighbor matching.


texts:
I've seen several questions asking how to convert .tab files to .shp files in QGIS, MapInfo etc, which is what I'd like to do. However I have several files and I was wondering if there was a tool available to do a batch convert of several files in a folder from .tab to .shp?

I've had a look around but not found an answer, apologies if this is a duplicate, if so feel free to just point me in the right direction.


texts:
I have a pretty large polygon shapefile (almost 1.5 million buildings) that I want to display in Google Maps. As the KML would be far too big, I have to create tiles, right?

I don't want to host my own server, but just upload the tiles on my site. But how do I create the tiles for such a shapefile?

TileMill seems not working with such a large shapefile. It has to be Google Maps as I want to use the geolocation service that works best for the country the data is from.

Any hints?


texts: Selecting annotation by color to correct Annotation Class using ArcPy?

texts:
I have some municipal data that was converted into a GIS and all of the address numbers were coded as dimensions in the annotation feature class. There is an annotation class for address, but none of the addresses are of that type. However, all of the addresses are displayed on the maps in a particular color (Mars Red/255 0 0), which makes me think there is a way to select these address numbers.

I have found that there is an arcpy tool called Select Features by Symbol Color (SelectFeaturesBySymbolColor_production) but I do not have the production mapping extension available.

Is there a way to make this selection using the mapping module, perhaps? Or an open source solution? Or if I could make a definition query for those it would be easy to quickly select them.

I'm using ArcMap 10.6


texts: How big areas/.osm data files can Maperitive handle?

texts:
How do the file-size limitations of Maperitive work out?

Background: I'm interested in making a map of the whole of Germany. Any pre-generated maps I find online only cover a region at a time, so maybe most mapping software has a limit.

Would it help if I make a .osm file with some features removed? Basically I only need streets and street numbers of civilized areas - so Schwarzwald might be ok to skip :-)


texts:
I am able to manipulate coordinates, their position and diverse data. when testing I can see the street name on the map but i can't seem to find the data property to get.

What i am trying to do is simply get the street name based on my feature coordinates.

I logged the feature but i am can't find the data unless i missed it. I checked the API too but there is nothing mentioning streets properties.

UPDATE

As suggested I used Nominatim reverse like this

https://nominatim.openstreetmap.org/reverse.php?format=json&lat=LAT&lon=LON&zoom=16

and through a simple xhttp request like this I got the response.

       var xhttp = new XMLHttpRequest();
    xhttp.onreadystatechange = function() {
        if (this.readyState == 4 && this.status == 200) {
            document.getElementById("demo").innerHTML =
            this.responseText;
       }
    };
    xhttp.open("GET", "ajax_info.txt", true);
    xhttp.send();

Source : W3C

Note I am using v3.11 and OSM.


texts:
I downloaded the data from here
https://geo-nsdi.er.usgs.gov/metadata/map-mf/2323/metadata.faq.html
When I loaded it into QGIS it says it has no CRS set. There is no .prj file.
After trying a variety of different CRSs I tried making my own custom CRS
The string I came up with is.
+proj=lcc +lon_0=-117.3819444 +lat_0=23.0 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +a=6378206.4 +lat_1=33.0
This is my first time trying to make a custom projection.
The data is put in Canada north of where it should be.
It should be in the Great Basin in the USA mostly in Nevada.
Here is the pertinent information I used to guess at the string.
Any ideas on what I can change to have it placed in the correct place?
What coordinate system is used to represent geographic features?
The map projection used is Lambert Conformal Conic.
Projection parameters:
Standard_Parallel: 33.0
Standard_Parallel: 45.0
Longitude_of_Central_Meridian: -118.0
Latitude_of_Projection_Origin: 23.0
False_Easting: 0
False_Northing: 0
Planar coordinates are encoded using Coordinate pair
Abscissae (x-coordinates) are specified to the nearest 130.0
Ordinates (y-coordinates) are specified to the nearest 130.0
Planar coordinates are specified in meters
The horizontal datum used is North American Datum of 1927.
The ellipsoid used is Clarke 1866.
The semi-major axis of the ellipsoid used is 6378206.4.
The flattening of the ellipsoid used is 1/294.98.
Vertical_Coordinate_System_Definition:
Altitude_System_Definition:
Altitude_Datum_Name: National Geodetic Vertical Datum of 1929
Altitude_Resolution: 30.0
Altitude_Distance_Units: meters
Altitude_Encoding_Method: Attribute values


texts: Installing httplib2 through the Python console in QGIS 2.18.0 on Windows OS

texts:
I want to install httplib2 through the Python console in QGIS 2.18.0 on Windows OS, in order to use another plug-in instagram2qgis.
I downloaded and unzipped the folder of httplib2 but I don't know how to install it.


texts: How to generate a point distance matrix using shortest network distance?

texts:
My problem is as follows.

To start with, I have the following spatial datasets:

A road network;
13,000 odd point locations representing populated localities (ie places where people live);
Service (point) locations (2000+) representing the location of a specific human service;

For each populated locality (origins), I need to find the network distance to the 5 nearest service locations (destinations), and record those 5 distance measures against the originating populated locality.

Software resources available to me include ArcMap (ArcView licence only), Manifold (with Business Tools extension), Quantum GIS, GvSIG.

Any help or suggestions welcome.


texts: Speed of Layer Rendering and Visibility in ArcGIS Online

texts:
This map is a combination of several layers: http://arcg.is/2fTGANA

When I share this out with others they experience difficulty in rendering when they log in to AGOL as compared to when they do not.  What could be causing this?

Also, even though I adjust the visibility on the raster layer it does not necessarily display correctly.


texts:
I am trying to use the spatial autocorrelation function in the statistical analysis toolbox but I am unsure what I need to set as the inputs. I want to use it to test for autocorrelation in a set of 'samples' from the seafloor that have been used to test the accuracy of an automated classification of substrate type.

The Input Feature Class would presumably be the locations of the samples as X, Y points in the map, and then is the Input Field the 'value' for the substrate class assigned to that point?

Another issue is that the 'samples' we are using are along transects, how accurate can an autocorrelation test be when the points are located along a line?


texts: How to spiderfy stacked POIs in CARTO's standard map?

texts:
I'm creating this map on CARTO.

Some POIs are stacked even at max zoom level, so some POIs are impossible to be clicked.

Are there any methods to make these POIs spiderfied?


texts:
A recent post on Cartastrophe refers to a "Worldwide Campaign to Eliminate Needless North Arrows."

So when, if ever, are North arrows necessary?

As an interesting aside the concept of North being up is not an arbitary one, but in the Northern hemisphere

  "due to the visible rotation of the
  night sky around the visible celestial
  pole"

however

  "Chinese and Islamic culture even
  considered south as the proper top end
  for maps."

http://en.wikipedia.org/wiki/North


texts: Does gwc support tile expiration configuration for arcgisLayer?

texts:
I am setting up a GeoWebCache standalone server, and  configed an ArcGIS layer with this configuration below .
<arcgisLayer>
    <name>default</name>
    <tilingScheme>C:\gwc\data_dir\data\default\conf.xml</tilingScheme>
    <tileCachePath>C:\gwc\data_dir\data\default\_alllayers</tileCachePath>
    <hexZoom>false</hexZoom>
    <expireCacheList>
       <expirationRule minZoom="0"  expiration="-2" />
    </expireCacheList>
    <expireClientsList>
       <expirationRule minZoom="0" expiration="7200"/>
    </expireClientsList>
</arcgisLayer>

According to the document, expireCacheList enable server-side "tile caching", expireClientsList enable client-side "tile caching" by setting http response header with right values, am I right?
But those getMap request got responsed with no-cache headers, why? How can enable client-side caching for arcgisLayer?
Cache-Contorl:no-cache Expires:0


texts: Distance between a line and a specified waypoint in lat long

texts:
I'm based in the UK and I'm trying to determine the distance in meters between a line (created via a set of lat lng waypoints) and a specific point (also defined in lat lng).I've had a look at a couple of the posts here and come up with the following code.The distance I'm looking for is the shortest distance between the line and the specified waypoint.

I'm using org.geottols and com.vividsolutions.jts implementation.

Can someone please verify if I've done this right?

    double distance = 0.0;GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();

    Point startPoint = geometryFactory
            .createPoint(new Coordinate(prow.getStartLatitude(), prow.getStartLongitude()));
    Point endPoint = geometryFactory
            .createPoint(new Coordinate(prow.getFinishLatitude(), prow.getFinishLongitude()));

    List<WayPoint> waypoints = getWaypoints(prow);

    int length = waypoints.size() + 2;
    Coordinate[] coords = new Coordinate[length];
    coords[0] = startPoint.getCoordinate();

    Point[] waypointList = new Point[length];
    waypointList[0] = startPoint;
    waypointList[length - 1] = endPoint;

    int i = 1;
    for (Iterator<WayPoint> iterator = waypoints.iterator(); iterator.hasNext();) {
        WayPoint wayPoint = (WayPoint) iterator.next();
        Point wp = geometryFactory.createPoint(new Coordinate(wayPoint.getLatitude(), wayPoint.getLongitude()));
        waypointList[i] = wp;
        coords[i] = wp.getCoordinate();
        i = i + 1;

    }

    coords[length - 1] = endPoint.getCoordinate();

    Point targetWaypoint = geometryFactory.createPoint(new Coordinate(waypointLat, waypointLng));

    LineString route = geometryFactory.createLineString(coords);
    route = (LineString) Densifier.densify(route, .01);

    try {

        CRSAuthorityFactory crsFac = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", null);
        CoordinateReferenceSystem wgs84crs = crsFac.createCoordinateReferenceSystem("4326");
        MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, wgs84crs);
        Geometry prowLine = JTS.transform(route, transform);
        Geometry waypoint = JTS.transform(targetWaypoint, transform);

        Double degrees = prowLine.distance(waypoint);
        // g3.isWithinDistance(g4, distance);

        distance = degrees  * (3.14159/180) * 6371;

        Coordinate[] c = DistanceOp.nearestPoints(waypoint, prowLine);

        Coordinate c1 = new Coordinate();
        // System.out.println(c[1].distance(g4.getCoordinate()));
        JTS.transform(c[1], c1, transform.inverse());
        // System.out.println(geometryFactory.createPoint(c1));
        distance = JTS.orthodromicDistance(targetWaypoint.getCoordinate(), c1, DefaultGeographicCRS.WGS84);

    } catch (Exception e) {
        e.printStackTrace();
    }

    return distance;


texts:
I'm doing a spatial query, but I always get an error:

  "msPOSTGISLayerGetShape(): Query
  error. msPOSTGISLayerGetShape called
  on unopened layer (layerinfo = NULL)"

This is my code:

layerObj layerUtbs = Util.MSMap.getLayerByName("Utbs");
if (layerUtbs.queryByRect(Util.MSMap, rect) == (int)MS_RETURN_VALUE.MS_SUCCESS)
{
    resultCacheObj resultado = layerUtbs.getResults();

    if (resultado.numresults > 0)
    {
        // here i got a error.
        layerUtbs.getFeature(resultado.getResult(0).shapeindex, resultado.getResult(0).tileindex);
    }
}

Here's my map file configuration

LAYER
        CONNECTIONTYPE postgis
        NAME "Utbs"
        CONNECTION "user=pmc dbname=SiteInovacao host=localhost password=123 port=5432"
        PROCESSING "CLOSE_CONNECTION=DEFER"
        DATA "geom from (SELECT *, oid FROM pmc.utbs) as utbs using unique oid using srid=31983"

...

Thanks!


texts: Subset a SpatialPolygonsDataFrame by ID in R?

texts:
I've used the gIntersect function from the rgeos package to determine which polygons intersect in two SpatialPolygonsDataFrames (spdfs). The result is a logical matrix with rownames comprising the polygon IDs from sfdf #1 and colnames the polygon IDs from spdf #2. I'd like to use this matrix to subset the two original spdfs, omitting those polygons that have no overlap, before using the function gIntersection. I can easily convert the logical matrix into two vectors of IDs (one for each spdf) that represent those polygons with at least some overlap.

My question is - can I use a vector of polygon IDs to subset a spdf? If so, will it correctly subset values in the @data slot? I know I can generally subset spdfs using this syntax:

spdf[spdf@data$myVariable[someVector], ]

This won't work for IDs, as they are buried much deeper inside the spdf list structure.


texts: What is the recommended workflow for using TileMill and OSM?

texts:
Using osm2pgsql, I imported an OSM extract into a Postgres database. Then I made a layer in in TileMill that contains the data from Postgres. The map was exported to a MBTiles file. At this point, all my transportation features (lines) are on a single layer and my map is terribly unpleasant to look at.

How do I apply the openstreetmap-carto cartocss to my map? I discovered I can create a style for each layer but that seems like a painfully slow way to style my map and I'd need to import each feature onto its own layer resulting in a separate query for each feature line. There's gotta be a better way.

What is the recommended workflow for building a map in TileMill from OSM data in a postgres db?


texts:
I have a polygon and I want to make parallel lines inside of the polygon. I need to have the maximum length when I sum the length of all of them.
Each line must be 15 m apart.


texts: calculate flow accumulation in a river with two outlet points

texts: I'm trying to calculate flow accumulation layer (with arcGIS) in a basin with 2 outlet points. I have used "dem reconditioning" tool to modify my dem according to a vector layer, which contains my digitized stream network. The problem is that a few km before the main river delta a diversion is constructed (for flood management) creating thus 2 riverbeds, the natural and the constructed one, resulting in two outlet points and two deltaic parts along the coastline. I'm trying to work with a specific flood event in the river and I need to calculate flow accumulation layer with these two outlet points . I've tried many things with no success.


texts: Avoiding multiple elif statements in Python parser for ArcGIS label expressions?

texts:
My code for labels look like this:

 def FindLabel ( [NAME_GEO] ):
  if [NAME_GEO]=="Tbilisi" :
     return "<BOL><FNT size ='14'>" + [NAME_GEO] + "</FNT></BOL>"
  elif [NAME_GEO]=="soxumi" :
     return "<BOL><FNT size ='12'>" + [NAME_GEO] + "</FNT></BOL>"
  elif [NAME_GEO]=="baTumi" :
     return "<BOL><FNT size ='12'>" + [NAME_GEO] + "</FNT></BOL>"
  else:
     return [NAME_GEO]

It works but as you can see there are two 'elif' statements for cities, which are pretty much the same. If I wish to have more city names with same label properties as it is for "soxumi" and "baTumi", I would need to add more 'elif' statements.
I wanted to change code, so there would be only one 'elif' statement for "soxumi", "baTumi" and any other cities, however I failed... could anyone suggest how to do it?


texts: Custom CRS settings disappears after file reloading

texts:
I have create custom CRS for my project. That's one of local CRS for Russia - MSK-72 / zone 1. Parameters are: +proj=tmerc +lat_0=0 +lon_0=66.08333333 +k=1 +x_0=1500000 +y_0=-6000000 +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs.
After switching project CRS to this one custom, all layers on my map are located on correct places. I've also reprojected layers to custom CRS. After reloading the project, layers still in correct position, but CRS both for project and for all layers is indicated as UNKNOWN. How can I fix this error?


texts: Snapping points to line yields error 000800?

texts:
I am trying to snap points (temp2.shp) to nearest line (Rail2.shp) at any of its vertex within any distance.

ArcPy code:

temp2 = "M:\\temp2.shp"
arcpy.Snap_edit(temp2, ["M:\\Rail2", "VERTEX", '10000 Feet'])

Error Message:

Failed to execute. Parameters are not valid.
ERROR 000800: The value is not a member of END | VERTEX | EDGE.
ERROR 000800: The value is not a member of END | VERTEX | EDGE.
ERROR 000800: The value is not a member of END | VERTEX | EDGE.
Failed to execute (Snap).


texts:
I am trying to figure out why pgr_trsp() isn't giving me the interpolated route (start and end positions are the full line instead of somewhere along that line).

I believe it has to do with the way I a joining on the edge table. My theory is that when I join on the edge table it uses the geometry stored in this table which doesn't have substrings of the line and so the whole line is selected.

Below is a screenshot of what I am trying to get.

With the pgRouting plugin option trsp(edge), it provides the red line as shown but I am getting the blue line.

Code below:

SELECT st_linemerge(edgeTable.geom_way) FROM pgr_trsp('SELECT id, source, target, cost FROM edgeTable',
(SELECT id FROM origin),
.5,
(SELECT id FROM destination),
.5,
false, false) AS shortestPath
JOIN edgeTable ON shortestPath.id2 = edgeTable.id;

Can anyone tell me how to join this properly so I can use the results from the function?


texts: Qgis gdalogr:contour is not generating output

texts:
I am trying to use gdal:contour script in Qgis console and then script. Here is code which I tried write to console.

raster = QgsRasterLayer('C:/Users/Bulva.user/Desktop/test_raster.asc', "DEM")
QgsMapLayerRegistry.instance().addMapLayers([raster])
processing.runalg("gdalogr:contour", raster, 1.0, "ELEV", None, "C:/Users/Bulva.user/Desktop/contour.shp")

the output of console is:

{'OUTPUT_VECTOR': 'C:/Users/Bulva.user/Desktop/countour.shp'}

but there is no shapefile as output. I tried processing.runandload and there is a problem with loading layer, because shapefile wasn't generated.

I tried it few days ago and I swear it worked but today I have no idea why it is not working.

This my log in Log Messages Panel > Processing:

2016-05-11T11:45:35 2   Problem with OTB installation: OTB was not found or is not correctly installed
2016-05-11T11:47:25 0   GDAL execution console output
            Usage: gdal_contour [-b <band>] [-a <attribute_name>] [-3d] [-inodata]

                                [-snodata n] [-f <formatname>] [-i <interval>]

                                [-f <formatname>] [[-dsco NAME=VALUE] ...] [[-lco NAME=VALUE] ...]

                                [-off <offset>] [-fl <level> <level>...]

                                [-nln <outlayername>] [-q]

                                <src_filename> <dst_filename>

            FAILURE: Too many command options.

and in OGR:

2016-05-11T11:47:25 1   Data source is invalid ()
2016-05-11T11:47:25 1   Data source is invalid ()


texts: Increasing speed limit by 35 percent in ArcGIS Field Calculator?

texts:
I have a network with a speed limit attributes.

I want to increase the speed limit by 35 percent in order to run an analysis for an ambulance.

For example, if the speed limit is 60, I will increase it by 35 percent.

Can I do this in the field calculator?


texts: Using CityGML data with ArcGIS Desktop or QGIS?

texts:
I'm looking for advice on how best to use the CityGML data-structure and GIS (preferably ArcGIS/QGIS). I've read through their website and wiki and saw a list of software here that can be used to view the data, but I don't know which one works best with GIS.

I have not used CityGML before (as I'm still figuring out what software to use to view the data) so I would appreciate any suggestions.

My goal is to use the CityGML approach (with varying levels of detail, or LODs) to perform spatial analysis at different levels of accuracy.


texts:
I want to connect Mapnik to MSSQL Spatial SDE, but I don't know how to do this.

  import mapnik
m = mapnik.Map(256,256)
mapnik.load_map(m, "test.xml")
m.zoom_all()
mapnik.render_to_file(m, "2.png")

xml code is:

 <Map background-color="#ffffff00" srs="+init=epsg:3857">
  <Style name="My Style">
    <Rule>
      <PolygonSymbolizer fill="blue" />
      <LineSymbolizer stroke="red" stroke-width="0.1" />
    </Rule>
  </Style>
  <Layer name="world" srs="+init=epsg:3857">
    <StyleName>My Style</StyleName>
    <Datasource>
<Parameter name="type">ogr</Parameter>
            <Parameter name="string">MSSQL:server=xxxx;database=xxxx;UID=xxxx;PWD=xxxxx</Parameter>
            <Parameter name="layer_by_sql"> SELECT top 10 Shape FROM dbo.xxxx</Parameter>
            <Parameter name="extent">0 ,0,0, 0</Parameter>

</Datasource>
</Layer>
</Map>


texts:
I installed PostgreSQL 9.6, PostGIS 2.3.1 and Pgadmin III on my computer, and I was reading the documentation, so I thought that would be a great idea see the function on practice and then modifie them to learn a little bit, but the most SQL codes that I see on documentation doesn't work if I put them on pgadmin3, for example on the page 130 in postgis's documentation (http://postgis.net/stuff/postgis-2.3.pdf) there is the following statement:

SELECT ST_MLineFromText(’MULTILINESTRING((1 2, 3 4), (4 5, 6 7))’);

When I put it on pgadmin3, I receive an error that says:

Error of sintax in or next to "2".

LINE 1: SELECT ST_MLineFromText(’MULTILINESTRING((1 2, 3 4), (4 5, 6...

                                                     ^

********** Error **********

ERROR: error sintax in or next to "2"
SQL state: 42601
Character: 45

By the way: I'm using Debian Sid

Sorry if it's a dumb question, I'm a newbie in spatial databases.


texts:
Running the code through Console Phyton in QGIS Desktop:
copy_renderer_layer = QgsProject.instance().mapLayersByName("map_layer_linestring")[0]
paste_renderer_layer = QgsVectorLayer("LineString?crs=EPSG:4326", "out_layer", "memory")
pr = paste_renderer_layer.dataProvider()
pr.addAttributes(copy_renderer_layer.fields())
paste_renderer_layer.updateFields()

paste_renderer_layer.setRenderer(copy_renderer_layer.renderer())

The code works correctly, but after executing it three times, the QGIS Desktop program simply closes! What is the problem?


texts: Debugging error 000865 in Python script tool for use with raster elevations?

texts:
I'm trying to create a python script tool that will take in a DEM/raster as a parameter, pull out the maximum and minimum elevation values, and iterate through this range of values at a specified counting interval.

This script outputs the resulting raster masks to an output geodatabase, as another parameter.  I keep getting the error '000865: Input Raster: inDEM does not exist' when I try to run my script tool.  I'd like help resolving this issue in my code, as I'm not sure what the cause is.

import arcpy
from arcpy.sa import *
import os

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

inDEM = arcpy.GetParameterAsText(0)
intervalElevation = int(arcpy.GetParameterAsText(1))
outGDB = arcpy.GetParameterAsText(2)

elevMinResult = arcpy.GetRasterProperties_management("inDEM", "MINIMUM")
elevMaxResult = arcpy.GetRasterProperties_management("inDEM", "MAXIMUM")

minElevation = elevMinResult.getOutput(0)
maxElevation = elevMaxResult.getOutput(0)

myElevations = range(minElevation, maxElevation, intervalElevation)

for elev in myElevations:

    outDEM = os.path.join(outGDB, "lessthan" + str(elev))
    outLessThan = LessThan(inDEM, elev)
    outLessThan.save(outDEM)
    print "Finished " + outDEM


texts: Importing a simple elevation raster to ArcGIS Online

texts:
I have a very simple elevation raster that I want to bring in to AGOL for use in a web map or web scene. I need to preserve the elevation data. I have tried to bring it in as a tile service but I lose all the elevation data. I feel like this should be a fairly simple task but I've spent hours trying to figure it out.
Any suggestions?


texts: I would like to know a way to measure performance of SOC machines. We are in process of moving from physical servers to VM Servers.


texts: arcPy SelectLayerByAttribute_Managment SQL Expression issues using wildcard

texts:
I have a Layer of Streams called "My_Channels" with the following source

Data Type:  File Geodatabase Feature Class
Database:   C:\Users\user\Documents\ArcGIS\Default.gdb
Feature Class:  My_Channels
Feature Type:   Simple
Geometry Type:  Line
Coordinates have Z values:  No
Coordinates have measures:  No

I Then Copy the layer so I don't destroy it and create "My_Test_Channels"

I'm attempting to use the arcpy.SelectLayerByAttribute_management({layer},{selection},{WHERE})

I'm executing the code populated with the following SQL WHERE CLAUSE

arcpy.SelectLayerByAttribute_management("MY_Test_Channels","NEW_SELECTION","UNIT_NO= 'E%'")

So far I've tried looking at the arcPy Docs for this tool.

I then want to copy the specific selections to a new layer and put that in my Default.gdb file. To be Clear I'm trying to extract Streams starting with the UNIT_NO "E". I have been able to execute the next code which is arcpy.CopyFeatures_management("My_Test_Channels","cropedLayer") The reason I would like to do these types of operations is I'd like to study those specific stream types and currently our layer is setup with all of the streams for our area."E" as their UNIT_NO There are other stream groupings with other letter designations

Am I using the wrong tool to do this since the feature is a line and not necessarily a layer?

Or is it simply that my SQL expression is not formatted correctly?

I've also looked at using Select Analysis Tool but am unsure which is easiest, I plan to do this on other datasets what may be polygons, polylines, lines, points, etc. The current license that I have does not include the Split tool.


texts: Opening dbf without losing ObjectID and geometry

texts:
I found a lot of related questions, but not completely an answer to my struggle. I can open the dbf-file in ACL and save it as dbf. I added some attributes and wanted to open it again in ArcGIS, but apparently I lose 2 columns when opening my dbf in ACL (and it does too when I open it in Excel): the objectid and the geometry.

How can I prevent this from happening or are there other options I have to join the object ID's with a new attribute that I only have in a table format (no GIS-coordinates)?


texts: Visualizing large number of WKT files in QGIS

texts:
I have a directory which stores the result of breaking up a big geometry into many parts. There are 100 files named part-0000.wkt to part-0099.wkt. I need to point QGIS to this directory and have it load and render all the wkt files. How can I do this? I read this, but the QuickWKT plugin can't even load geometry from a file, much less a directory.

Also I am trying the add delimited layer tool and the OK button never gets enabled! What's up with that?


texts: How to use google maps engine connector?

texts:
I am trying to use google imagery in QGIS through the Google maps engine connector plugin.  I relatively new to QGIS and need some help.  What I want to do i use the google maps in QGIS and then print them on a custom 3ft by 3ft map for a client.  Is this possible with the plug in?  Also if it is how do I get it to work.  I have established an API as directed but now I am stuck.

Is there anymore step by step documentation on how to use the plugin?


texts: What's the best way to add Tilemaps to QGIS without having to use Webmercator?

texts:
I'm looking for a good way to add tilemaps in qgis that are either in a custom projection or that I can reproject in qgis.  I have used the openlayers plugin to add basemap layers, but it forces me to use webmercator.  I have also spent some time playing with adding layers as a gdal wms as described here http://www.gdal.org/frmt_wms.html, but layers added in this manner are very slow even in their native projection.

My end goal would be to have ortho photos in a fast format in my local projection, and to be able to make some custom basemaps.


texts: How to Generate Points on Corner of Polygon Only

texts:
I have a polygon layer made up of grid cells, illustrating a survey system:

I'm trying to create a point on each corner of the polygon. I've used a geoprocessing tool called "Feature Vertices to Points" in ArcGIS, but what it does is create a point for each vertice:

Is there a tool out there that will allow for the creation of corner points?

Please note this is a small subsection of the survey grid that I clipped out to run a proof of concept.

What I'm basically trying to do is create these points of this grid (which was created many years ago) and then apply the same methodology to the latest version of the survey grid, which is based off of NAD83, to see how much and where the grids have shifted.

The differences in the grids are important when determining how measurements change when it comes to locating buildings, wells etc.


texts: How to Hide Toolbars from QGIS1.7 programmatically?

texts:
I am a beginner to QGIS 1.7 and getting problem in hiding the toolbar.
I am using:

fileTB = self.iface.fileToolBar()
fileTB.setVisible(False)

This does not work but doesn't give an error either.

I also tried:

self.iface.mainWindow().removeToolBar(fileTB)

This permanently removes the toolbar

fileTB.clear()

but I don't want that, because I want it to be get invisible, when the user click the 'on' button it should get displayed.

Update:

def initGui(self): # I used this one
    fileTB = self.iface.fileToolBar()
    fileTB.setVisible(False) # As You suggested me qgis.utils.iface.fileToolBar().setVisible(False)


texts: Obtaining Non-Summary Zonal Statistics from Raster Data

texts: I am using ArcMap 10.5.1. The data I am looking for are zonal statistics from a raster using polygon-defined zones. I can use the Zonal tools to calculate summary statistics, but I am looking for the entire range of values. Surely, if it can calculate the summary stats, it must be able to create an output of some format with the entire range of individual values. Any ideas? (FYI: I am unacquainted with writing code for ArcGIS)


texts: Does anyone know of a good unaffiliated review of Rockworks?

texts: I've been struggling to find a good software package for drillhole data management, subsurface geology modelling, and visualization (both 2D and 3D): something that can also be used to produce sophisticated vertical and inclined cross sections.  I took the Gemcom GEMS course and found the software too complex (database setup, plotting) and rigid. It's also astronomically expensive. I evaluated Target for ArcGIS by Geosoft using their tutorial and data and found it lacking features, buggy and slow. (Dumbed down version that works with ArcGIS.) As I was to evaluate the full standalone version of Target I stumbled over a recommendation for Rockworks by RockWare. The videos on their site look amazing. Too good to be true to be honest. I would like to read some unaffiliated reviews of the software but am having a difficult time finding anything valuable. Can someone please provide a link a share a few thoughts about the software package?


texts: Why do GDAL and rasterio give different projected bounding coordinates for a given input geotiff?

texts:
I have a geotiff, and I want to find its extents and reproject them. I can use rasterio:

# Destination CRS
dest_crs = some CRS

# Get projected boundaries with rasterio
clipbox = rasterio.open('mygeotif.tif')
cl, cb, cr, ct = clipbox.bounds
l, b, r, t = transform_bounds(clipbox.crs, dest_crs, clipbox.bounds[0],
clipbox.bounds[1], clipbox.bounds[2], clipbox.bounds[3])
print(l,b,r,t)

-75.0277857229692 -9.469504586908355 -74.09777685862646 -7.796173999420965

...or I can use GDAL:

# Get projected boundaries with GDAL
clipbox = gdal.Open('mygeotif.tif')
wktproj = clipbox.GetProjection()
sref = osr.SpatialReference()
sref.ImportFromWkt(wktproj)

# Get its extents
ulx, xres, xskew, uly, yskew, yres  = clipbox.GetGeoTransform()
lrx = ulx + (clipbox.RasterXSize * xres)
lry = uly + (clipbox.RasterYSize * yres)

# Reproject the extents
transform = osr.CoordinateTransformation(dest_crf, sref)
tl, tt, garbage = transform.TransformPoint(ulx, uly)
tr, tb, garbage = transform.TransformPoint(lrx, lry)
print(tl, tb, tr, tt)

-75.02766365052823 -9.468344577976415 -74.09777685862646 -7.797125769469675

So one latitude value is the same, but the other three are different. Does anyone have any ideas why this would be the case? I should note that cl == ulx, cb == lry, etc. I.e. the input boundary coordinates are the same pre-projection.

EDIT:

Pre-projection coordinates:

print(ulx,uly,lrx,lry) # GDAL

496950.0 -861870.0 599040.0 -1046760.0

print(cl, ct, cr, cb) # rasterio

496950.0 -861870.0 599040.0 -1046760.0

Projections:

print(dest_crs) # GDAL

"osgeo.osr.SpatialReference; proxy of Swig Object of type 'OSRSpatialReferenceShadow *' at 0x117b85e10"

and using a function called wkt2epsg, we see

print(wkt2epsg(wktproj))

EPSG:32618

print(dest_crs) # rasterio

EPSG:32618


texts: How to convert ESRI .prj WKT string to Proj4 string WITHOUT GDAL

texts:
ESRI shapefiles include a .prj file containing CRS information in Well-Known Text (WKT) format.
I'm developing a software in which I need to systematically read shapefiles' CRS information, convert it to a Proj4 string, and feed Proj4 with it.

The software I'm developing can not have dependencies on GDAL.

So, how can I systematically convert/parse a WKT string from a .prj file into a Proj4 string? Is there any direct mapping between the tags?

Alternatively, instead of converting to a Proj4 string, finding the corresponding EPSG code would also be acceptable.


texts: Plugin processAlgorithm. Same code, different results in script vs plugin code

texts:
I'm adding functionality to a QGIS Plugin to create a print layout and add stuff to it like a map, legend, etc.

To test this code, I wrote a script. The script works and creates the print layout the way I want. Now I'm trying to add the code into the processAlgorithm() of my plugin's code. I'm getting different results with the same code from the script.

Here is what I have so far:

    def processAlgorithm(self, parameters, context, feedback):
        """This actually does the processing for creating the print layout and exporting as .pdf"""
        #needs all the arguments (self, parameters, context, feedback)

        log = feedback.setProgressText

        input_vector = self.parameterAsVectorLayer(parameters, self.INPUT_VECTOR, context)
        input_title = self.parameterAsString(parameters, self.INPUT_TITLE, context)

        log(f"Input vector name: {input_vector.name()}")
        log(f"Title: {input_title}")

        """This creates a new print layout"""
        project = QgsProject.instance()
        manager = project.layoutManager()
        layout = QgsPrintLayout(project)
        layoutName = "PrintLayout"

        layout = QgsPrintLayout(project)
        layout.initializeDefaults()
        layout.setName(layoutName)
        manager.addLayout(layout)

        """This adds a map item to the Print Layout"""
        map = QgsLayoutItemMap(layout)
        map.setRect(20, 20, 20, 20)
        #Set Extent
        canvas = iface.mapCanvas()
        map.setExtent(canvas.extent())                  #sets map extent to current map canvas
        layout.addLayoutItem(map)
        #Move & Resize
        map.attemptMove(QgsLayoutPoint(5, 27, QgsUnitTypes.LayoutMillimeters))
        map.attemptResize(QgsLayoutSize(239, 178, QgsUnitTypes.LayoutMillimeters))

        results = {}
        return results

When I run this part of the code in the script, I get a nice map in my print layout. It looks like this:

When I run my plugin with the code above, I get this:

This map is the correct size and in the correct position. But obviously there is no map and something is broken.

Here is the log after I run the plugin:

Is the problem in the way I return the "results" from the processAlgorithm? In the QGIS documentation it says this:

However, I don't want to return anything. I just want to create the print layout with items on it. There is no output from this part of the plugin.

I just create "results" variable as an empty dictionary and return it because I am supposed to return a dictionary. So obviously I am missing something.


texts:
My Arcgis Desktop 10.1 was working fine until, suddenly, I became unable to use the Arctoolbox tools :

when I try to open them from the arctoolbox window they just do not answer,
when I try to open them from the search tool I get a warning/error message:

  Warning

  Adding this item is not supported or the application that opens it is not found.

I've done researches on internet, but although some problems deal with arctoolbox dysfunctionalities, I found none that are quite like mine.

Anybody has a suggestion or solution?


texts: Spatial statistics within a polygon in Google Earth Engine

texts:
I am wondering if the information provided in this exchange is still accurate.

number of pixels inside a polygon in Google Earth Engine

I am trying to calculate mean precipitation within a river basin, where precipitation data are in a raster and the basins are a polygon. Using Google Earth Engine, I have already calculated average areal weighted precipitation, but I am having trouble discovering the exact formula/calculation process GEE uses to calculate an areal weighted mean.

According to the GEE Changelog on 2017-08-28, there appears to be no minimal coverage required for a pixel to be added to a weighted calculation. Is this correct? So, even if 0.01% of a precipitation pixel is contained within a basin, then it should be included in the calculation of a weighted mean?


texts: Areas of polygon within different buffer overlaps using ArcMap

texts:
I have two shapefiles, one with some 100 buffers and one with polygons of which I want to calculate the areas within the buffer overlap areas, but with regards to the adjacent buffers. E.g. the area sum of the polygons within the overlap between 2 and 4, 3 and 2 etc.
So far, I have calculated the overlap zones by using the intersections tool and then intersect those areas again with the polygon shapefile. Summing up manually is a very tedious job, but I cannot get my head around how to achieve this automatically with ArcMap.
Edit: I maintained all the original attributes of  both shapefiles (centroid, unique ID etc.) Each buffer has an unique ID, which has been transfered to the clusters as well. I have duplicates within the intersection shapefile depending on the number of buffers which intersect. Same applies to the polygons whose areas I want to calculate.


texts: Creating empty feature class in file geodatabase based on schema of ArcSDE feature class?

texts:
I have 150+ feature classes in SDE for which I need to create an empty feature class in a file geodatabase. For these feature classes I need to import the attributes and the name of the feature class as in SDE. I tried to export to an XML workspace document but it is not working because of some error on the SDE and I'm not an SDE administrator, so I am unable to modify the things on the SDE.

An alternative to this that I tried was creating a script in Python for creating the feature class. I get the names of the feature class from the SDE but am unable to create them on the local SDE.

Also how can I read the geometry type for the feature class in arcpy?

Here is the code I've written:

import arcpy
arcpy.env.workspace= r"C:\Users\088927\AppData\Roaming\ESRI\Desktop10.0\ArcCatalog\ENT-QA @ sample.sde"
fcList = arcpy.ListFeatureClasses('sde."ENT-QA".*',"","")
output = "C:/Users/088927/Desktop/Schema.gdb"
for fc in fcList:
    fullName = arcpy.ParseTableName(fc)
    #fcTemplate = fc.
    finallist = []
    nameList = fullName.split(",")
    databaseName=nameList[0]
    ownerName = nameList[1]
    fcName = nameList [2]
    fcStr = str(fcName)
    print fcStr
    arcpy.CreateFeatureclass_management(r"C:\Users\abc\Desktop\Schema.gdb",fcstr)


texts:
I'm trying to get the area of a raster (without NULL cells) using gdalinfo, I read the gdal documentation but found no parameter that would give me the area; I can calculate it from the histogram obtained, by subtracting the values in the first histogram bin (which where the NULL in this case) from the width*height values (size) and multiplying by the resolution, but hist function takes some time and maybe NULL cell will have another value in other raster files. This is the function.
gdalinfo -json  -hist raster.tif | jq ((.size[0]*.size[1]-.bands[0].histogram.buckets[0])*.geoTransform[1]*.geoTransform[1])/10000'


texts:
I have a geodatabase that I would like to convert to shapefile with epsg:4326 projection. For this I use ogr2ogr.

However when I reproject, the results are off. They seem to have shifted?

Any idea why?

The commands I'm using:

/opt/anaconda3/envs/python35/bin/ogr2ogr -f 'ESRI Shapefile' /volumes/data/Y2018M08D13_RH_Process_Basisregistratie_Gewaspercelen_V01/process2_V07/BRP_Gewaspercelen_2009.shp /volumes/data/Y2018M08D13_RH_Process_Basisregistratie_Gewaspercelen_V01/process_V07/BRP_Gewaspercelen_2009.gdb

and then convert to epsg 4326 (I used to do everything in one step but split for debugging purposes)

/opt/anaconda3/envs/python35/bin/ogr2ogr -t_srs EPSG:4326 /volumes/data/Y2018M08D13_RH_Process_Basisregistratie_Gewaspercelen_V01/output_V07/BRP_Gewaspercelen_2009.shp /volumes/data/Y2018M08D13_RH_Process_Basisregistratie_Gewaspercelen_V01/process2_V07/BRP_Gewaspercelen_2009.shp

Geodatabase: http://geodata.nationaalgeoregister.nl/brpgewaspercelen/extract/2009-definitief/brpgewaspercelen.zip

The file geodatabase has this metadata crs:

+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.87740,4.0725 +units=m +no_defs


texts: How to measure area of polygons generated from gps data in qgis?

texts:
Somebody classified my question as a duplicate. But applying the answer doesn't work for me. Here is what I did:
- I made sure that my layer is in a projected coordinate system (UTM30N, WGS84)
- I used the field calculator
- I used the area command.
- I get 0.
- And I don't know where to adjust the units...

I have gps data from around 200 farm areas and I want to compute the correct area of each polygon. What I did so far is this:

I imported the data from around 10 garmin devices to basecamp, gave them IDs and exported them in one gpx-file (around 2 MB). Than I imported the GPS data to qgis and used the "line to polygon" tool to generate shapefiles.

My coordinate system is WGS 84 and UTM 30N for the south west of Burkina Faso. The whole project in qgis as well as the layers have this coordinate system.

When I create a new column in the attribute table, open the field calculator and use the "area" command, I get a "0" for each polygon. When I use the measure tool I get very small figures (f.e. 0.024 m2 instead of 43.635 ha).

I use QGIS 1.7.0 (Wroclaw) on a MacBookPro Mac OS X (10.5.8).

Best, Martin


texts: Why is the OSM map itself is not visible in my qgis window?

texts:
I'm new in working with QGIS.

I have tried to download an OSM map of Rome and its surroundings into my QGIS. As far as I am concerned, I have followed each step which I could find on various websites. In fact I do see the three layers of my OSM map in the side bar (points, polylines and polygons). However, I do not see the map itself.

Does anyone have an idea what my mistake might be?

In the meantime, I tried another way. I have downloaded the "Isle of Man"-Shapefile from Geofabrik. I was able to import all vectors for all data that are meant to appear in the map. In the end I saw the outline of the Isle of Man, the houses, the rivers etc. However, I do not see the map itself.

I have a new laptop, 64-bit, Windows 8, QGIS 2.4


texts: Editing outline color when drawing quantities using graduated colors in ArcMap?

texts:
I have constructed a basic map of population density for a given group over an area, and I am creating a keyframe to show the overall area with my extent indicator labelled.

The problem is that the grey outlines between city zones at such a large scale become visual clutter, and I'd like to remove them from that layer. I'm having difficulty finding a way to alter the color (and select no color) for the outlines in that keyframe layer since I'm using graduated colors, where the option to alter outline color/width is not available.

Does anyone know a way to do so using ArcMap 10.1?

Here's the web of grey borders


texts: Converting multiple lines into single line using ArcGIS Desktop?

texts: I've no idea if this is basic or what but I have this data that contains multiple lines (forming a polygon shape). I want to convert those lines (multiple - each with a start and end point) to a single line. I tried dissolve, multipart to single part etc. I cannot use topology since it influences precision. I have zero knowledge of programming.


texts: Excel (*.xls) imported into File Geodatabase is missing records?

texts:
I am having a few issues with a join to a File GDB from an Excel .xls file I imported as a table to the FGDB.

The problem being records that are in the original spreadsheet are not appearing in the imported table and when I join and keep all records I am getting null values although I know in the original spreadsheet the values exist. There about 65000 records in the spreadsheet.

I have also tried with a CSV and failed as well as changing integer and double in the field type - can't seem to work out why this is.

The join is one to one, so nothing complex about it at all.

I'm using ArcGIS 10.2.2 for Desktop.


texts: Intersection error QGIS - invalid geometries that are actually valid

texts:
I'm trying to intersect two vector layers (having the same CRS EPSG:4326 WGS 84, in QGIS 2.18.10) but I get this error in the log:

Input layer A contains invalid geometries (feature 52).
Unable to complete intersection algorithm.

I have checked the validity of both layers and in both cases there are no errors or invalid geometries. Nevertheless, I have deleted feature 52 of the layer with the invalid geometry, and when I intersect them again, I find an error on feature 51, and so on until I delete every single feature.

Is there a way to check for errors in the geometry other than 'check validity'?
Could it be another type of error that QGIS is trying to indicate?


texts:
So there is the cassini projection in Proj. It is possible to use it in the OSGeo4WShell.
example: +proj=cass +lon_0=.... +lat_0=.... +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs
Is there a way to define such CRS in Python?
example for ecef system: ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
But I want to use the cassini projection (to transform to ecef).
How do I define it?
Is CAssini implemented in pyproj?


texts: Does pgRouting support turn-by-turn routing?

texts:
I have route list as result of pgRouting with street names. How to get messages "Turn left", "Go along" etc ...
Can pgRouting do this or I need some other source and which?


texts:
I test cartodb with some data and I found it very useful. I found this instructions for installing cartodb locally in github What is CartoDB? and I try to understand if cartodb is a CMS for maps.

Is like the joomla but for maps?

I mean can I install it and use it like the joomla?

The https://cartodb.com/ web site is just a host service for maps?


texts:
I am working with  Landsat 5 images in UTM-WGS84 using the Qgis 2.14 software.
I need to obtain the XY pixel coordinates with their ND and to extract them in an Excel list.
How I can do it?


texts: Merging features using Collector for ArcGIS (Field Maps)

texts:
In 'Collector for ArcGis',  is it possible to merge two polygon features (from two different layers) into one feature?
It would be required to decide from which layer the attributes are used.


texts: How to retrieve GIS data in Multiple Orthophoto at point level

texts:
I have been using a python script that can read Data from Orthophoto(*.tiff) of size 300Mb-2Gb. But this is taking 32s per request of lat-long, which is not ok.

I have different Orthophotos with Different data like DSM, DTM, NDVI etc of same location at different times.

What would be an ideal way or tools to do this? This has to be done as a program, as the this needs to be replicated for various locations.

Currently I am evaluating storing this data in DB against lat-long, with all values of DSM, DTM, NDVI etc across time being associated with the lat-long key in mongo or PostGIS.

Is there a better way of solving this issue?


texts: Measuring distance between discharge point, along river, until next intersection with another river?

texts:
i hope someone can answer me my simple question. I have two sets of data

a river network --> line feature
discharge points along those rivers --> Point Data

The discharge points represent wastewater treatment plants. I try to find those treatment plants which don't have a inflow of another river, for a minimum of 2 km downstream. This means the river from this point on does not intersect with another river for at least 2km.

How can i find those plants? Are there any tools for this? What about the Arc Hydro tools?


texts:
I frequently am in need of color gradients that aren't readily available through ArcGIS or the ColorBrewer add-in. I find the process of manually creating gradients in Arc very tedious: saved custom colors don't come up in the edit color ramp window.

I am wishing that there was some way to use Python to batch generate colors ramps for ArcGIS. For example, if I gave it a "high" color and a "low" color, it could create and import a linear gradient between the colors. Importing might simply be a matter of saving to a specified directory.


texts: Unable to implement pgrouting(plugin) pgr_astar routing algorithm (shortest path)

texts:
I have been working on a project for indoor navigation.

I made the floorplan of the building and road-map for routing purposes as shp files (using JOSM, making GeoJSON and then converting to shp file using http://www.mapshaper.org/)

Then using POST GIS ShapeFile Loader (of pgAdmin III) I imported road-map(route.shp) to the database.

In QGIS, i was able to use the shortest path plugin to do the routing

//////////////////////////////////////////////////////////////////////////////////////////

EDITED

//////////////////////////////////////////////////////////////////////////////////////////
The following are screen shots of route

The following are screen shots of route_vertices_pgr

When I tried to use the pgrouting plugin of QGIS , I got the below error

What am I doing wrong?


texts: How to use the ExactIntersect flag for QgsFeatureRequest() object

texts:
I have two layers, a polygon_lyr and a line_lyr. I want to use a QgsFeatureRequest() for the line_lyr.getFeatures() but instead of returning the lines which intersect the bounding box of a polygon, I just want the lines which intersect the actual polygon geometry. I assume this is what ExactIntersect flag does. But I do not know how to actually use it.

for i, poly in enumerate(polygon_lyr.getFeatures()):
    poly_geom = poly.geometry()

    for line in line_lyr.getFeatures(QgsFeatureRequest(poly_geom.boundingBox())):
        #processing operations


texts: If I've given two different layers the same CRS, why won't they overlay one another?

texts:
I have a "towns" and "state" shapefile for Massachusetts as two separate vector layers in QGIS. I've set the project coordinate system to WGS 84 at the bottom of the screen (next to the "coordinate" and "scale" windows), I've also used the "define current projection" tool under Vector> Data Management Tools. Lastly, I've right clicked on each shapefile in my layers window to "Set Layer CRS." I've also enabled "on the fly" transformations. Is there a common reason why the two shapefiles will not sit atop one another, though I've given them the same projection?

Thanks!
Alex


texts: QGIS displays satellite data very blurry when ArcGIS Pro looks clear?

texts:
I am looking at a WorldView 3 image in QGIS (right in image) and in ArcGIS Pro (left) and there is a huge difference in quality. QGIS looks so blurry/a lower resolution.

Any idea what's happening?


texts: What is the correct way to query for signboards related access tags in Overpass in an area?

texts:
I am new to all this stuff.
How can I search (i.e., on Overpass) which access tag can be added to the segments on OSM, based on sign board information I have. Is there a way to search access tag related information on overpass based on sign boards?


texts:
I've written a python script to plot a cross-section from a DEM along a line.

Everything works fine except when I Show the plot, ArcGIS hangs and I have to force it to close with Task Manager.

I am working with Windows 7 and ArcGIS 10.3. Following is my code snippet:

import arcpy
from arcpy import env
import matplotlib.pyplot as plt
...
...
...
plt.plot(x,y)
plt.xlabel('Distance (m)')
plt.ylabel('Elevation (m)')
plt.title('X-section along the line')
plt.grid(True)
plt.savefig('Xsection.png')
plt.show()
plt.close()


texts: real world coordinate when digitizing raster image to vector in QGIS

texts:
I'd like to create a vector layer by digitizing from an image.  How can I get the (nearly) correct longitude and latitude when digitizing from an image from Google Map and from an image from a satellite?

For example, I would like to digitize the following image.  I know a lat- lng pair of one location and I know the scale.  How can I create the vector layer to match this?  What CRS should I use?


texts: High values when opening HDF bands in rasterio

texts:
Im getting very high values when opening a band in rasterio from an HDF file. I got the same thing when I open netcdf files. I guess the issue for netcdf is related here: (1) and (2)
Thus when opening a netcdf, I apply scaling and offsetting.
ds = rasterio.open(f'netcdf:H08_20210101_0000_1MARP031_FLDK.02401_02401.nc:AOT_L2_Mean')
band = ds.read(1) * ds.scales + ds.offsets
band = band.astype('float32')

Should I do the same for HDF file since it's the same issue?
ds = rasterio.open(
    f'HDF4_EOS:EOS_GRID:"MOD11C3.A2018001.006.2018032091151.hdf":MODIS_MONTHLY_0.05DEG_CMG_LST:LST_Day_CMG')
band = ds.read(1) * ds.scales + ds.offsets
band = band.astype('float32')


texts: Slope raster seems to come out distorted?

texts:
I've been trying to get a slope raster from a DEM, so I can reclassify it and aggregate it with land use data.

The problem is that my slope raster seems to have a lot of distortion (diagonal lines). I've tried changing some settings but I really have no idea what can be causing this. I'm using a projected CRS, EPSG 31983. I've got the DEM from whole state where I live in Brazil, then I had to clip it to my AOI, and that is all the processing I made with the DEM before trying to get the slope.

Some information about the DEM file from the metadata:

  Driver: GDAL provider
  GTiff GeoTIFF

  Dataset Description
  /Users/Guilherme/Desktop/FINAL/DEM/DEM_SIRGAS_CLIPPED.tif

  AREA_OR_POINT=Area
  DataType=Generic

  HighCellSize=0.039888
  LowCellSize=0.000277
  MaxCellSize=1.994400

  Band 1
  SourceBandIndex=0

  Dimensions - X: 692 Y: 346 Bands: 1

  Origin: 227546,7.41239e+06

  Pixel Size: 28.8724,-28.8724

  No Data Value: 0

  Data Type: Float32 - Thirty two bit floating point

  Pyramid overviews

  Layer Spatial Reference System: +proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

  Layer Extent (layer original source projection) -
  227546.4710596329532564,7402400.2644091211259365 : 247526.1508292471116874,7412390.1042939284816384

  Band No: 1
  Min Val: 531.9765014648
  Max Val: 960.2605590820
  Range: 428.2840576172
  Mean: 608.2492439474
  Sum of squares: 1395234469.1412367821
  Standard Deviation: 76.3367047628
  Sum of all cells: 145634332.9768066406
  Cell Count: 239432

DEM file:

Slope generated from DEM (a lot of odd diagonal lines):

EDIT:

I went to the source again, and got the raw DEM to process/clip it again. This time around everything seemed to work just fine. I believe what happened was some kind of mismanage when warping the raster.

I didn't get the chance to try Freighty's solution because SAGA and Orfeo are not working properly for me, but I believe if you have access to them, it's a feasible solution.


texts: Calculating Euclidean distance in ArcGIS

texts: I am willing to model the distribution of an amphibian species in relation to climate and land cover in MaxEnt. For this I have downloaded the WorldClim data and the Global Consensus Land Cover (https://www.earthenv.org/landcover). The latter data consists of a raster file for each land cover class (open water, cultivated and managed vegetation,...etc.). However, I will not be using directly the raster of each land cover class, but instead a raster of the Euclidean distance of the XY coordiantes (longitude, latitude) where the species occurs from the nearest open water, cultivated and managed vegetation,...etc. I tried to calculate the Euclidean distance in ArcGIS 10.7.1 with the following: input raster (shapefile deature class of my occurrence data) and input barrier raster (open water raster). The Euclidean distance is calculated without errors but I get a blank raster with only one value (!) and I can not find where the problem is. Any thoughts?


texts: Batchprocess asking for CRS for every file

texts:
I "wrote" a model using the graphical modeller which uses GDALs "translate" to convert one geotiff to one georeferenced jpeg. The model simply tells, which is the desired CRS and - most important - omits all otherwise manually necessary settings.
I then start the model as a batch process, since the conversion has to be done 350 times, loading all the to be converted files. For me it works fine.
For a colleague on the other hand, who would like to use this model, QGis asks for every file for a CRS.

Do you have any suggestions, why this could happen and how we could fix it?


texts: Automatically save outputs variables in folder workspace using ArcPy?

texts:
I want to create python script with aim to import him in environment arcmap 10.1 and to make a hydrologic analysis

because the scipt meybe can have many outputs raster file or feature class would like to find a way for automatic save outputs variables in the folder workspace

for example

I want to import a raster file (Dem) and to take as export a 3 raster file from the command fill sinks, Fill Sinks,flow directions and flow accumulation from arc hydro toolbox

cl = my import dem
fill = output for Fill Sinks
fdr = output for flow directions
Fac = output for flow accumulation

and i want export them(3 raster) automatically saved to the folder workspace without definitive output parameters from the user

import arcpy

arcpy.ImportToolbox("C:\Program Files (x86)\ArcGIS\Desktop10.1\ArcToolbox\Toolboxes\Arc Hydro Tools.tbx", "archydrotools")

# Local variables:
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = arcpy.GetParameterAsText(0)
cl = arcpy.GetParameterAsText(1)
# Process: Fill Sinks
arcpy.FillSinks_archydrotools(cl, fill, "", "", "ISSINK_YES")
# Process: Flow Direction
arcpy.FlowDirection_archydrotools(fill, Fdr, "")
# Process: Flow Accumulation
arcpy.FlowAccumulation_archydrotools(Fdr, Fac)

any idea ?


texts: R language, problems with SpatialPixelsDataFrame

texts:
(This question comes from https://stackoverflow.com/questions/40698369/r-language-problems-with-spatialpixelsdataframe?noredirect=1#comment68629240_40698369)

The following two scripts will generate a "SpatialPixelDataFrame" object:

# FIRST
library(rgdal)
elev.grid <- readGDAL("whatever.asc")
elev.grid <- as(elev.grid, "SpatialPixelsDataFrame")

# SECOND
library(raster)
library(SDMTools)
library(adehabitat)
elev.grid <- raster("whatever.asc")
elev.grid.asc <- asc.from.raster(elev.grid)
elev.grid.SPDF <- asc2spixdf(elev.grid.asc)

However, the first one exceeds the capability of my computing resources when applying it to big (15000 x 16000) grids, and the second one generates an object which I can't use for some of my further analyses. For example, when I use it for krige purposes

x <- krige(V3~var, points, elev.grid)

I get the following:

  Error in model.frame.default(terms(formula), as(data, "data.frame"), na.action = na.fail) : invalid type (closure) for variable 'var'

Whether providing me a trick to bypass the memory/capability issue in the first case (preferably), or fixing the error generated by the second case please recommend.