texts: How to create polylines from a point Feature Class linking all points using arcpy

texts:
I want to create polylines linking a bunch of points. Individual polylines must be created between every point. I started by reading the x & y cords and then created unique pairs as well as all combinations. Unfortunately my python is lacking and I'm also stuck with ArcGIS 10.0.

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

myt = ()
newt = ()
pairt = ()

# Get the XY coords from the point featureclass
inputfc = 'C:\\Temp\\testing.gdb\\BunchofPoints'
desc = arcpy.Describe(infc)
shapename = desc.ShapeFieldName
rows = arcpy.SearchCursor(inputfc)
for row in rows:
    feat = row.getValue(shapename)
    pnt = feat.getPart()
    myt = myt +((pnt.X),(pnt.Y))

#Create pairs
newt =zip(myt[::2],myt[1::2])

#Create combinations
for pair in itertools.combinations(newt,2):
    pairt = pairt + pair

#Check
print pairt

#Create lines
mypoint = arcpy.Point()
myarray = arcpy.Array()
featList = []
for feature in pairt:
    for coPair in feature:
        mypoint.X = coPair[0]
        mypoint.Y = coPair[1]
        myarray.add(mypoint)
    polyline = arcpy.Polyline(myarray)
    myarray.removeAll()
    featList.append(polyline)
arcpy.CopyFeatures_management(featureList, 'C:\\Temp\\all_lines.shp')

I get the following error:
    Runtime error : 'float' object is unsubscriptable

I'm struggling with how to read the pair of coordinates I've created and also with creating the individual lines.


texts: The next big thing in standardized spatial data exchange format?

texts:
I remember reading something a few months back, in the latter half of 2012, about a coordinated group effort to design a new standard for sharing geographic / spatial data. It has or aims to have the potential of replacing the showing-its-age (20+ yrs) de-facto standard of the Esri shapefile as the means of choice for widely sharing our data across all devices and platforms. Aside from addressing the most pressing shapefile limitations, this new thing incorporates raster and metadata too.

This big-important-thing I don't recall is not SpatiaLite - Shapefile of the future? or Are there any attempts to replace the shapefile? but some of the conversation swirling around it was similar. I spent quite some time today trying to track this down, but my search-fu proved not up to the task (which might not bode well for its future).

Please, what is the name of this elusive project that has the potential to deeply impact all our workflows? (Assuming it gets off the drawing board, lord knows we need something to.) And, how or where does one get involved or at least keep track of it?


texts:
I'm referring to this icon

I'm just curious as to its origins and what the name of the symbol is. If memory serves, Google Maps introduced it and the other mapping libraries adopted it. Would anybody know what the symbol is called or is it just a "map marker"? Did it exist before web maps or is purely a product of the web map age? Or is just a stylized pin?


texts:
I have a simple Geodatabase in access format. I'd really like to host this database on an ArcGIS Server through one of their services.

The idea is that my users would be able to connect to this Geodatabase to look at both their features, and their databases, and then relate them. I do not want to allow editing of the database, just viewing purposes.

If at all possible, I'd really like to stay away from using ArcGIS Engine Runtime on the client side.

What might be the best service I can run? And does anyone have a good example of how to connect to that service to extract the database information?


texts:
I'm having some trouble getting this map to work in IE8. Currently tried using some of the old IE8 CSS for both but that didn't seem to work. Any ideas for this?

Here is the code for the main map:

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/themes/css/cartodb.css" />
<link rel="stylesheet" href="css/style.css" />
<link href="css/leaflet-search.css" rel="stylesheet" type="text/css" />
<link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css" />
<!--[if lte IE 8]>
<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v2/themes/css/cartodb.ie.css" />
<![endif]-->
<!--[if lte IE 8]><link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.6.4/leaflet.ie.css" /><![endif]-->
<script src="js/cartodb.js"></script>
<script src="js/leaflet-search.js"></script>
<script type="text/javascript" src="js/leaflet-hash.js"></script>
<script>
var map;

function init() {
  // initiate leaflet map
  map = new L.Map('map', {
        zoomControl: true,
        center: [49.2500, -123.1000],
        zoom: 10,
  })
   L.tileLayer('https://a.tiles.mapbox.com/v3/cardinalmaps.jf059aad/{z}/{x}/{y}.png', {
    attribution: ''
  }).addTo(map);

  var layerUrl = '';
  var sublayers = [];
  var LayerActions = {
    reset: function(){
      sublayers[0].setSQL("SELECT * FROM do_not_update_bc_shelter_data");
      return true;
    },
    yearround: function(){
      sublayers[0].setSQL("SELECT * FROM do_not_update_bc_shelter_data WHERE open_closed = 'Year-Round Shelter'");
      return true;
    },
        closed: function(){
      sublayers[0].setSQL("SELECT * FROM do_not_update_bc_shelter_data WHERE open_closed = 'EWR Shelter - Closed'");
      return true;
    },
    openshelters: function(){
    sublayers[0].setSQL("SELECT * FROM do_not_update_bc_shelter_data  WHERE open_closed = 'EWR Shelter - Open'");
      return true;
    },
  };

  cartodb.createLayer(map, layerUrl)
    .addTo(map)
    .on('done', function(layer) {
      // change the query for the first layer
      var subLayerOptions = {
        sql: "SELECT * FROM do_not_update_bc_shelter_data",
      }
      var sublayer = layer.getSubLayer(0);
      sublayer.set(subLayerOptions);
      sublayer.on('featureClick', function(e, latlng, pos, data, layer) {
        console.log(pos)
        map.panTo(latlng, {animate:true})
      });

      sublayers.push(sublayer);
    }).on('error', function() {
      //log the error
    });

  $('.button').click(function() {
    $('.button').removeClass('selected');
    $(this).addClass('selected');
    LayerActions[$(this).attr('id')]();
  });

      map.addControl( new L.Control.Search({
        url: 'http://nominatim.openstreetmap.org/search?format=json&q={s}&countrycodes=ca&state=bc&bounded=1',
        jsonpParam: 'json_callback',
        propertyName: 'display_name',
        propertyLoc: ['lat','lon'],
    zoom: 13,
      }) );
}

And here is the map online:

http://cardinalmaps.com/files/bchousing-4/


texts:
I've got Collection 2, level 2 Landsat 8 dataset which contain mountains. The problem is the cast shadows from the mountains, so that I can't classify these pixels well. These pixels have much lower values because of the shadows.
Is there any method to remove these shadows?

I'm using ArcMap, and I've already tried the Topographic Correction, but it wasn't successful.


texts: Usage of "160" instead of "156" as the UN M49 or ISO 3166-1 numeric country code for China

texts:
I'm guessing the short answer is "Something complicated and political to do with Taiwan", but I'd like to know the detail...

I'm matching up some different international data sources and ran into a small snag because China didn't match up. One data source (some health data published by the UN agency IARC) had China as "160" (and had no entry for Taiwan in their data set) while every other source I've seen, including the official ISO list and UN M49 list, has China as "156" (with Taiwan as "158", and "160" undefined).

For every other country, this IARC data matches the "ISO 3166-1 numeric" or "M49 UN" standard, China seems to be the only anomaly. I can't find anything explaining this alternate standard.

The only clue I can find is their "Data Sources And Methods" page for China (country=160 in the URL...), which includes a map that depicts Taiwan as part of China.

Is 160 a lesser-used country code for "China, and Taiwan"? What's its official status?

I've tried looking on the official ISO 3166-1 pages but couldn't find anything relevant.


texts: I'm using SQL Server 2008 R2 to store Geometry data types.  I have a bunch of lines as individual features, that share common names ('NAME' field).  I would like to create a view that Groups By 'NAME' and combines the geometries.  I can only find STUnion function and I'm not sure how to use this to in a query to do what I want.


texts: Contour Plugin Processing Algorithm is not detected in PyQGIS 3 Standalone

texts:
I was trying to access processing algorithm using PyQGIS 3 on a standalone script. I found this solution and tried it but when I list the processing algorithm using print(QgsApplication.processingRegistry().algorithms()) it only lists the built-in plugins. The contour plugin that I wanted to import is not listed although it works just fine QGIS.

I have tried this code in QGIS Console:

import processing
QgsApplication.processingRegistry().algorithms()

and the contour plugin is listed there as <contour.ContourGeneratorAlgorithm.ContourGeneratorAlgorithm object at 0x7fe052e6aca8>, but when I tried in standalone script using this script:

sys.path.append('/usr/share/qgis/python/plugins/processing')

import processing
from processing.core.Processing import Processing
Processing.initialize()

print(QgsApplication.processingRegistry().algorithms())

it only lists the built-in plugins. Where is my mistake?

Btw, Contour plugin is an external QGIS plugin that make a contour from data points as an input. This plugin still in experimental state.


texts: QGIS not exporting JSON to GeoJSON format

texts:
My PostgreSQL/PostGIS table have a JSON column, how to export it as GeoJSON properties to QGIS?
That is, I have a PostgreSQL connection at QGIS and it is reading all other tables, except this one that have a JSON colummn. It is  possible to read JSON column at QGIS?
NOTES
Community say that "QGIS was born for PostgGIS"... Seems a bug because JSON  good support exist since PostgreSQL v9.3, and GeoJSON is used by QGIS since 2016 or before, it is not a strange object...
If it is not a bug, and it is not an expected QGIS feature, a workaround is perhaps to create a SQL VIEW where each line is a GeoJSON to be exported to QGIS... I am testing  "generic VIEW template" solution, using jsonb_populate_record function.


texts: Clipping imagery to regions of interest with Google Earth Engine before exporting

texts:
My code currently collects all Landsat imagery and then creates a mean summer NDVI for each year.  I would like to clip the the summer raster to an vector polygon shapefile (such as attached).

Currently, I think my code just filters full Landsat rasters that have at least partial overlap with any polygon element in the ROI leaving me with huge, multigigabye rasters to export.  I only want the raster pixels within my ROI.
I would like to clip my combined Landsat 4-8 collections after line 120.
https://code.earthengine.google.com/620d62244bfe1cf910d7c38ddcd65187
edit:
Here is an asset for designating my ROI.
https://code.earthengine.google.com/?asset=users/edeutsch/ROI_2016


texts: Passing two variables into codeblock of Calculate Field tool of ArcGIS Desktop?

texts:
I have a script tool that truncates each word in field to 3 characters and outputs the result in a new field:

"THIS IS AN EXAMPLE STRING" becomes "THI IS AN EXA STR"

import arcpy
arcpy.env.overwriteOutput = True

# Prelogic used to pass through for calculating truncated name values
codeblock = """
def nmshort(stname):
  s = stname.split()
  s2 = s[:] = [elem[:3] for elem in s]
  return ' '.join(s2)
""" # Finished code block

try:
    # Set input parameters
    inFC = arcpy.GetParameterAsText(0)
    inName = arcpy.GetParameterAsText(1)
    inNameShort = arcpy.GetParameterAsText(2)

    # set expression for field calculation to variables
    expression = "nmshort(!{0}!)".format(inName)

    # perform field calculation
    arcpy.AddMessage("Calculating truncated strings for values in input name field")
    arcpy.CalculateField_management(inFC, inNameShort, expression, "PYTHON_9.3", codeblock)
    arcpy.AddMessage("Calculation complete")

except:
    arcpy.AddError("Could not complete the calculation")

    arcpy.AddMessage(arcpy.GetMessages())

It works fine except I want to make a parameter in the tool that will allow the user to specify the length for truncation instead of a constant of 3.  With having to pass the codeblock to calculatefield as a string, I'm struggling to figure out how to pass the length parameter to a variable in the codeblock.  Most of the time the length will be either 3 or 4, so I suppose I could create two blocks of code: one to run if the length parameter is 3 and another if 4.  That would be the easy way out though and I want to learn (I'm a noob).  I'm really interested in doing this with an update cursor, if it can be done that way.  I've played with cursors with some success with simpler calculations.  I can't figure out how to write the more advanced logic involved in this calculation though.

I'm using ArcGIS 10.0 so I can't use the data access cursors.


texts: mixing lon/lat and mercator layers in a MapServer .map file

texts:
I have a basic question about my MapServer .map file, which has two layers, each with a different coordinate system, longlat vs. mercator.

See my .map file below. If I leave it as shown below, with the "version 2" UNITS/EXTENT active, it correctly draws my "MrSid" in the right place on the earth,  but doesn't draw my "shape layer" layer. On the other hand, if I comment out the "version 2" UNITS/EXTENT and activate the "version 1" UNITS/EXTENT, it then draws "shape layer" but not "MrSid".

So it's acting like the the UNITS/EXTENT has to match those of the layer for the layer to be drawn. Is that right? If so, how do I mix layers with different UNITS?

(By the way, I also have a secondary question about this .map file. If I don't include any PROJECTION entry at the outer/MAP level, it will draw the "shape layer" layer but not the "MrSid" layer. But if I add a PROJECTION, regardless of whether the projection is "proj=longlat" "ellps=WGS84" "datum=WGS84" or "proj=merc" "lon_0=0" "k=1" "x_0=0" "y_0=0" "a=6378137" "b=6378137" "towgs84=0,0,0,0,0,0,0" "units=m", it will draw either layer -- but not both at the same time, as per the main issue above. So my secondary question is why it seems to be happy using a default projection with my longlat shape file, but not with my mercator MrSid file.)

MAP
  IMAGETYPE   PNG

  # version 1
  # UNITS dd
  # EXTENT      -180 -90 180 90

  # version 2
  UNITS meters
  EXTENT -20037508.34 -20037508.34 20037508.34 20037508.34

  SIZE        550 450
  INTERLACE OFF
  PROJECTION
    # "proj=longlat" "ellps=WGS84" "datum=WGS84"
    "proj=merc" "lon_0=0" "k=1" "x_0=0" "y_0=0" "a=6378137" "b=6378137" "towgs84=0,0,0,0,0,0,0" "units=m"
  END
  SHAPEPATH   "../data"

  OUTPUTFORMAT
    NAME png
    DRIVER "GD/PNG"
    MIMETYPE "image/png"
    IMAGEMODE RGB
    FORMATOPTION  INTERLACE=OFF
    TRANSPARENT OFF
  END

  # Start of LAYER DEFINITIONS-------------------------------
  LAYER
    NAME "shape layer"
    DATA dourokukan
    STATUS DEFAULT
    TYPE LINE

    CLASS
      NAME "“shape class"
      STYLE
        COLOR 227 227 127
      END
    END
  END
  LAYER
    NAME "MrSid"
    DATA "c:\mrsid\adams\ortho_1-1_1n_s_nd001_2010_1.sid"
    STATUS DEFAULT
    TRANSPARENCY 100
    TYPE RASTER
    PROJECTION
      "proj=utm" "ellps=GRS80" "datum=NAD83" "zone=13" "units=m"
    END
    LABELITEM ""
    TEMPLATE VOID
  END
  # End of LAYER DEFINITIONS -------------------------------
END


texts: Text imported into SQL Server using ogr2ogr had extra characters at the beginning

texts:
I've managed to import some MapInfo data sets into SQL Server using ogr2ogr

But when I run a query in SQL Management Studio (SSMS), the various text fields look to be "blank"

However if I use the SUBSTRING function to trim off the first couple of characters, e.g. SUBSTRING(NAME, 3, 20), the text is displayed when I run the SQL query


texts: Seeking Public API to get GPS satellite's position?

texts:
Is there any public availible API that allows you to retrieve the current coordinates of satellites such as GPS and Galileo?

I've only found ways to get Satellites Ephemeries and calculate the coordinates from that, but those positions are up to 3 hours old. To solve this you've got to use a satellite future position predictor but the only API for that I've found only allows a projection of 300 seconds, so to make predictions up to the current time you've got to wait potentially for 4-5 minutes.

I'm thinking of making my own API that has everything already loaded and calculated to give the current satellite's position. I just don't want to reinvent the wheel if someone else has already done this.


texts: Understandiung difference between TIN interpolation methods Linear and Clough-Toucher of QGIS

texts:
I started to study TIN in QGIS and created my first TIN from a point layer with elevation. There are two interpolation methods called Linear and Clough-Tocher that give very different results. Linear looks like a regular Raster and with Clough-Tocher the result cannot be interpreted at first glance. Further the values of both methods vary a lot, as can be seen in the picture.

I cannot find information in the QGIS Docs what the difference is between the two methods.
What is the difference and why are the results so different?
I am using QGIS 3.4.13.


texts: Making ellipse geometry from some known parameters using PostGIS?

texts:
How to create an ellipse geometry with postgis from known axis coordinates and peak radius ?

From picture below :

Point A and Point B is known lat/lon

R is a known result of fresnelZone calculation (in Meters).

Yes I read this and this (using shapely instead of Postgis), but I don't know how to implement/adopt it to my need.


texts: PRACH all counters plot not well visualized

texts:
Here's the output that Im getting:

Im not sure why it is doing this..

I want it to be more like this

When i try plotting only two counters (two radii):
It is giving me this:

I hope it is clearer for you... The code that im using is the same that i posted here:
PRACH data mapped right but arcs not on top properly

Also, if i want to sort the percentages by descending order... How can i do that?


texts:
What is multiline string in gis? I have found when I use pgrouting certain edges have this geometric property mentioned as multiline string.

I am not sure what's the purpose of this geometry and how it actually defines an edge. Any insights?


texts: Raster algebra (subtract) changes extent of output

texts:
I have a tool that calculates volume of a portion of a raster DTM. I first clip a portion of the entire raster then fill it. I then subtract the clipped original from the filled raster. My problem is that if my project only has the DTM as part of it, it works fine. However if I include a much larger DTM in the Table of Contents, it doesn't, even though the larger raster is not referenced by the code. What happens is that the difference raster now has the extents of the larger raster and the differences are now near the upper left corner, instead of being aligned with the clipped & filled rasters. So under the first condition, the difference raster is 462x365, but in the second case, it balloons to 39000x39000, the extent of the second DTM.
I have included a set extent statement, but that seems to make no difference. Below is the offending part of the code. OK up to rastercalc_temp, but that raster has the aberrant dimensions.
Input_DEM = arcpy.GetParameterAsText(0)
Aguada_rectangular_boundaries = arcpy.GetParameterAsText(1)
Aguada_geodatabase = arcpy.GetParameterAsText(2)

# Set Geoprocessing environments
tempEnvironment0 = arcpy.env.extent
arcpy.env.extent = Aguada_rectangular_boundaries
arcpy.env.scratchWorkspace = "C:\Users\William\Desktop\ModelBuilder\Scratch\Aguada_scratch.gdb"
arcpy.env.workspace = "C:\Users\William\Desktop\ModelBuilder\Scratch\Aguada_scratch.gdb"

# Local variables:
Clipped_raster = "Clipped_raster"

# Process: Clip
arcpy.Clip_management(Input_DEM, "", Clipped_raster, Aguada_rectangular_boundaries, "-3.402823e+038", "ClippingGeometry", "MAINTAIN_EXTENT")

# Process: Fill
Fill_NCALM = arcpy.sa.Fill(Clipped_raster)
Fill_NCALM.save("Filled_raster")

rastercalc_temp = Raster("Filled_raster") - Raster("Clipped_raster")
rastercalc_temp.save("rastercalc_temp")

I have noticed that Clipped_raster and Fill_raster differ by 1 pixel depending on the options for Clip_management. However, I played around with the various possibilities with the same result.


texts: I need to download some cities boundaries polygons and store them in a Postgis database. Is this possible? Using what tools?


texts: What is the fastest way to select oldest, unique geometry from a table containing duplicates in PostGIS?

texts:
I have an output from an ArcPy arcpy.Union_analysis operation and for the overlap polygons that were created I need to select out the oldest feature (along with the attributes).  The result is 75740 rows.

Normally, I would just use FindIdentical_management, but I found the results unreliable for 10.2.2 (and I can't change the version I'm on). After rolling my own in arcpy.da and not getting the performance I wanted, I am trying to select out via a PostGIS query, using this answer to Delete duplicate geometry in postgis tables as template.

While this is a big improvement, I am not as experienced with PostGIS + window functions and I think there might be a faster approach?

So far, I have:

-- dtime , numeric lower values are older
WITH unique_geoms (
    rowid
    ,geom
    )
AS (
    SELECT row_number() OVER (
            PARTITION BY geom ORDER BY dtime ASC
            ) AS rowid
        ,geom
        ,oid
        ,dtime
        ,dval
    FROM x_dupes
    )
SELECT rowid
    ,geom
    ,oid
    ,dtime
    ,dval
FROM unique_geoms
WHERE rowid=1;

Explain: http://explain.depesz.com/s/9cV

"CTE Scan on unique_geoms  (cost=10312.46..12016.61 rows=379 width=104) (actual time=701.221..2049.404 rows=40389 loops=1)"
"  Filter: (rowid = 1)"
"  Rows Removed by Filter: 35351"
"  CTE unique_geoms"
"    ->  WindowAgg  (cost=8797.66..10312.46 rows=75740 width=154) (actual time=701.200..1577.202 rows=75740 loops=1)"
"          ->  Sort  (cost=8797.66..8987.01 rows=75740 width=154) (actual time=701.171..970.353 rows=75740 loops=1)"
"                Sort Key: x_dupes.geom, x_dupes.dtime"
"                Sort Method: quicksort  Memory: 24671kB"
"                ->  Seq Scan on x_dupes  (cost=0.00..2659.40 rows=75740 width=154) (actual time=0.011..272.713 rows=75740 loops=1)"

http://explain.depesz.com/s/9cV

PostGIS Version: USE_GEOS=1 USE_PROJ=1 USE_STATS=1


texts:
I am trying to publish TIFF in GeoNode. The TIFF file has a .tfw file.

Uploading directly from GeoNode fails. Then I uploaded to server and published it in GeoServer and it works in GeoServer. I want the layer to be available in data panel in GeoNode. For this I did geonode updatelayers but it gives ('Failed to process geotiff_coverage', AttributeError("'Coverage' object has no attribute 'srid'",))

Then I saved the image as GeoTIFF using QGIS (i.e. without .tfw file) but still same problem.


texts: Why do certain fields import as strings in QGIS?

texts:
I have imported a CSV to QGIS with lat/lon values representing points.  The other fields, at least in the original CSV file, contain numeric values.

> sapply(pune_hotspots, class)
         id         lat         lon        type        name         830         930        1030
  "integer"   "numeric"   "numeric" "character" "character"   "numeric"   "numeric"   "numeric"
       1130
  "numeric"

Now I would like to change the size of the points proportional to values in the attribute table, however, something peculiar happens when I open the Styling page under Layer Properties > Style

When selecting Categorized, I am able to choose any field in my point layer.

Here, I see that a number of fields have imported as strings(ABC), but they are actually filled with continuous values (see below)

Therefore , when selecting Graduated, I am only left with 4 choices of fields to represent, however, these are not the fields that I need to use.

My question is, when importing the CSV, how can I avoid the conversion of certain fields from numeric to strings in QGIS?


texts: How to integrate Live layer effect in QGIS 2.10 with the help of python code based automatation process?

texts: I need to work with live layer effect (blur, outer glow, shadow) etc . I know how to use this with the help of QGIS symbol tool. But i need to integrate this in to my python automated map generation process. I am searching code for this but unable to find out. If anyone knows it then please let me know


texts: Add shp attribute data to column with python script

texts:
If feature column (feature is the name of the column) in table is ferry crossing in the attribute data then I need to add Yes to the Ferry column that I have already added to the table.

If the feature column does not say ferry crossing I need to add No to the ferry column. Here is a code snippet that I am having an issue with.

I'm new to python.

#need to update FERRY column with yes or no values if ferry crossing is available in feature field
delimfield = arcpy.AddFieldDelimiters(fc, "FERRY")
cursor = arcpy.da.UpdateCursor(fc, ["FERRY"])
for row in cursor:
    if row[0] == "ferry crossing":
        cursor.updateRow("FERRY")
        print "YES"
    if row[0] is not "ferry crossing":
        cursor.updateRow("FERRY")
        print "NO"

del row
del cursor

#if feature is ferry crossing == YES in FERRY field
#if feature is not ferry crossing == NO in FERRY field


texts: Query about How to display arrow on road in two ways if road is twoway

texts:
I am a beginner with GeoServer, i built a code for road mapping and it work but now i want to display arrow on road according to road direction in both way if road is bidirectional,then arrow should display in both way.....

I have tried following code but it display the arrow only at center.....
please reply me for both direction....

    <Rule>
      <ogc:Filter>
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>FUNC_CLASS</ogc:PropertyName>
          <ogc:Literal>1</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>
      <TextSymbolizer>
        <Label>
          <ogc:Literal>;</ogc:Literal>
        </Label>
        <Font>
          <CssParameter name="font-family">Lucida Sans</CssParameter>
          <CssParameter name="font-size">15</CssParameter>
        </Font>

        <LabelPlacement>
        <LinePlacement>
        <PerpendicularOffset>-5</PerpendicularOffset>

        </LinePlacement>
        </LabelPlacement>

        <Fill>
          <CssParameter name="fill">#585858</CssParameter>
        </Fill>

        <VendorOption name="followLine">true</VendorOption>
        <VendorOption name="forceLeftToRight">false</VendorOption>
        <VendorOption name="forceRightToLeft">true</VendorOption>

      </TextSymbolizer>
    </Rule>
     </FeatureTypeStyle>


texts: Using Symbology to show count of points in Buffer with ArcGIS Desktop?

texts:
I'm pretty new to GIS and I'm using ArcGIS 10.4

I've gathered point information about two different variables, X and Y. I have created the buffer zone around the X points and used a spatial join to join the Y points within a .25 mile buffer around the X points.

How can I show the count of Y points inside this buffer using graduated symbols preferably?


texts: Using cost path in ArcGIS Spatial Analyst?

texts:
I'm new to using the Least Cost Path tool.

I have a depth raster (tiff format) and I would like to calculate the cost path between two points close to the land. I need to "force" the tool to calculate the path below 25 m and I do not know how to do that.

I have attached a pic where I show the cost path but I want it close to the coast.

Is there any way to do that?


texts:
In QGIS, I represent a polygon of the El Niño regions on the Pacific Ocean, as seen below:

Now, what I need is that the regions of Oceania and East Asia are shown on the left side of the map, similar to this image:

The shapefile of continents does not repeat itself, so how could it be done?
The idea is to avoid the use of secondary images like the OpenLayers plugin.


texts: Openlayers-3 wfs-t save coordinates displacement?

texts:
I'm using openlayers-3 in client side and geoserver in server side. I have a vector layer. I get this layer from geoserver as follow:

var url = "http://localhost:8080/geoserver/wfs?&service=WFS&" +
           "version=1.1.0&request=GetFeature&typeName=points:cinema" +
           "&outputFormate=application/json";

var format = new ol.format.WFS();
source = new ol.source.Vector({
    url: 'proxy.cgi?url=' + encodeURIComponent(u),
    format: format
});

cinemaLyr = new ol.layer.Vector({
    title: 'Cynema',
    source: source
});
map.addLayer(cinemaLyr);

I get cinemaLyr correctly, I modify this layer and save it as follow:

var node = wfs.writeTransaction(null, [modifiedFeature], null, {
    featureType: "points:cinema",
    featureNS: "http://opengeo.org",
    srsName: "EPSG:4326"
});
var s = new XMLSerializer();
var url = "http://localhost:8080/geoserver/wfs";

Ext.Ajax.request({
    url: 'proxy.cgi?url=' + encodeURIComponent(url),
    method: 'POST',
    xmlData: s.serializeToString(node),
    headers: {
        'Content-Type': 'text/xml'
    },
    success: function(response) {
        console.log("Everything is fine");
    }
});

This code save the modifiedFeature on server, but its coordinate is displacement(x save as y in geoserver and vice versa).  In Insert and Delete it do correctly, Only in Update request displacement point coordinate. How do I fix it?


texts: Finding closest feature to given address in ArcGIS API for Python?

texts:
I'd like to find the closest open space (i.e. a park, square, playground) from a City of Boston MapServer using Python's ArcGIS package.

I'd like to use arcgis.features.use_proximity.find_nearest() method. To that end I've geocoded the address and passed it first to a Feature constructor and then added that Feature to a FeatureCollection.

My ArcGIS public account does not support spatial analysis.

Does anyone have any recommendations for finding the nearest open space?

get_open_spaces.py

from arcgis.gis import GIS
from arcgis.features import Feature, FeatureCollection, FeatureLayer, FeatureSet
from arcgis.features.use_proximity import find_nearest
from arcgis.geocoding import Geocoder
import logging, os, requests

def _get_open_spaces():
    """
    Gets the open spaces from the City's MapServer
    """
    open_spaces_url = "http://gis.cityofboston.gov/arcgis/rest/" \
        + "services/EnvironmentEnergy/OpenData/MapServer/7"
    f = FeatureLayer(url = open_spaces_url)
    return f

if __name__ == "__main__":
    username = ************
    passwrd = *************
    log_filename = "log/open_spaces_log"
    try:
        os.remove(log_filename)
    except OSError:
        pass
    logging.basicConfig(filename=log_filename, level=logging.DEBUG)
    my_address = "94 Sawyer Ave., Boston, MA, 02125"
    g = GIS(username = username, password=passwrd)
    logging.debug("Using address: %s" % my_address)
    logging.debug("getting open spaces feature layer")
    near_layer = _get_open_spaces()
    assert near_layer, "something went wrong with the near_layer"
    logging.debug("now we're trying to get the address into its on feature_layer")
    coded_address = Geocoder(my_address, g)
    address_as_feature = Feature('point', coded_address)
    analysis_layer = FeatureCollection(address_as_feature.as_dict)
    find_nearest(analysis_layer, near_layer)


texts: i have developed the application for editing the map features online but i want the features to be saved as version of the layer not in original layer.this work is done using GeoExt +GeoServer+OpenLayers.


texts: how can I know which projection has this data?

texts:
I have data in an unknown projection to me.

I open it in QGIS and I see that the shapes are ok, but I need to identify the projection in order to merge this data with other stuff.

The points look like WGS84 (EPSG:4326) but without comma.

If a point in WGS84 is like this: (-58,-34)

In this projection the coordinates are something like: (5670000,6120000)

Any ideas?

EDIT:
Here there are three pics, showing the situation.

This is the data I'm trying to import, that i don't know what is the projection. Look at the point marked, it corresponds to the central square.

This is my data (it is in lonlat even though it looks deformed). The coordinates shown are from the same central square.

This is google maps showing also the central square coordinates.

SOLUTION image, transformed from Campo Inchauspe EPSG:22195 to lonlat. Thanks to mkennedy


texts: How to identify the upstream sub basins contributing towards different drainage outlet points?

texts:
I have done the watershed delineation of a  basin and have created several sub-basins using ARC-hydro tool.
Now I need to identify which are the correponding upstream sub-basins contributing towards a different drainage outlet point and find the total area of all sub-basins contributing towards a particular drainage point.
Is there any tool to perform this other than visual inspection because I have to deal with  a very large basin


texts: Setting origin raster to resample using R?

texts:
I need to resample a raster from nrow= 467, ncol=805, res=0.0416667 degree to nrow=2335, ncol=4025, res=0.00833333 degree. Basically make the pixel size smaller, because I have too many small polygons for raster statistics.

I have tried it with the code below, the problem is that the resulting raster is slightly above the origin raster, see figure.

In gray "test raster", in yellow "original raster".

The new image, raster comparison.

library(raster)
library(sp)
library(rgdal)
library(ggplot2)
library(itertools)

path = '~/path'
spath = '~/path1'

source <- raster(file.path(path,'file.tif'))
>source

Summary of origin raster:
class      : RasterLayer
dimensions : 467, 805, 375935  (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -119, -85.45833, 13.54167, 33  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : ~/path/file.tif
names      : raster_name

>summary(source)

        raster_name
Min.            5.516960
1st Qu.         9.151464
Median         10.334464
3rd Qu.        12.149574
Max.           17.560150
NA's            0.000000

# How the raster should be
target <- raster(nrow=2335, ncol=4025, extent(source), crs(source))
#origin(target) <- c(-119, 32.9583) # Commented for further 'experiments'
>target

class      : RasterLayer
dimensions : 2335, 4025, 9398375  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -119, -85.45833, 13.54163, 32.99997  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

> summary(target)
        layer
Min.       NA
1st Qu.    NA
Median     NA
3rd Qu.    NA
Max.       NA
NA's       NA

# Project source raster to target raster by bilinear
dest <- projectRaster(source, target, method='bilinear', tolerance=0.0001)
>dest

class      : RasterLayer
dimensions : 2335, 4025, 9398375  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -119, -85.45833, 13.54163, 32.99997  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : raster_name
values     : 5.517167, 17.58864  (min, max)

> summary(dest)
        raster_name
Min.        5.517167e+00
1st Qu.     9.154074e+00
Median      1.034819e+01
3rd Qu.     1.217350e+01
Max.        1.758864e+01
NA's        5.232403e+06

# Write new raster file to a path
writeRaster(flip(dest, direction='y'), file.path(spath, 'target.tif'))```


texts: Is it possible to hatch format railroads in leaflet?

texts:
Within mapbox tilemill it is possible to make make rail-roads like this:

(source: mapbox.com)

To achieve this you need to write the following code: Source

#railroads {
  ::line, ::hatch { line-color: #777; }
  ::line { line-width:1; }
  ::hatch {
    line-width: 4;
    line-dasharray: 1, 24;
  }
}

Would it also be possible to achieve the same within leaflet?


texts:
How do I build a legend in ArcGIS Pro, working with a Global Scene view?
If I work in a normal map, then in the "Insert" tab there is a legend function, and if I work with the globe, then this tab is not there.


texts: OpenLayers 4 WFS add new features in real time

texts:
I would like to refresh the WFS source of my Vector layer in OpenLayers 4 automatically. I have linestrings is stored in PostGIS 2.3 (it has 2 field: ID(pk) and geom(geometry)). The data is updated in real time. New inserts and sometimes updates (ST_AddPoint) in about every half seconds or more often. The data is published with Geoserver 2.12.
I'm using OpenLayers 4 to show the data on map.

Now, the problem is that the OL reloads the data continuously when the WFSrefresh function is called, so the lines are "flashing".

The goal is to not refresh the whole layer, just add the newest inserted objects to the map after the SQL insert (or update). How can I do this feature in OpenLayers?

The code:

        var vectorSource = new ol.source.Vector({
            loader: function (extent) {
                $.ajax('http://localhost:8082/geoserver/test/ows', {
                    type: 'GET',
                    data: {
                        service: 'WFS',
                        version: '1.1.0',
                        request: 'GetFeature',
                        typename: 'wfs_geom',
                        srsname: 'EPSG:3857',
                        bbox: extent.join(',') + ',EPSG:3857'
                    }
                }).done(function (response) {
                                        vectorSource2.addFeatures(formatWFS.readFeatures(response))
                });
            },
            strategy: ol.loadingstrategy.bbox,
            extractStyles: true,
            projection: 'EPSG:3857'
        });
            function WFSrefresh() { vectorSource.clear(); vectorSource.addFeatures(vectorSource.readFeatures()) }


texts: How can I correct/alter/change the workspace path in my environment settings?

texts:
This has been bothering me for a very long time.  No matter what I do, whenever I open ArcGIS Desktop, the environment settings for current and scratch workspaces revert to a path that is non-existent on my computer.

The blurred portion of the path, should be my username, but it is someone else's username.  If I change this, it reverts back when I reopen ArcGIS.  I'm trying to determine if there is a file somewhere (perhaps an xml file) which stores these default settings.  I tried deleting the normal.mxt file thinking that was the culprit, but to no avail.

It is getting very frustrating because whenever I try to run a tool, it defaults to this location and I have to change it manually.

I'm using ArcGIS 10.3.1 and Windows 7, 64-bit.

I'd prefer not reinstalling Arc, but I'm open to it if that is the only way.


texts: Why are the special fonts in a .mxd-project not used/found?

texts:
I have a .mxd-project file from a customer which uses special fonts (TrueType). I installed these fonts in Windows but my application (mostly a MapControl) doesn't use these fonts for displaying. I checked the mxd file (see whuber did here) in a text editor and the font names do appear.

I didn't find anything about using uncommon fonts in projects. I assume that it should work with simply installing the fonts in the OS but maybe there is something different to do?

What do I have to do to get the fonts displayed in my MapControl?


texts: GeoServer 2.15: "REST PathMapper Root directory path" ignored?

texts:
I am trying to create stores in my GeoServer setup through the REST API and try to force the store creation to take place in a specific directory by setting the "REST PathMapper Root directory path" under Settings > Global > REST PathMapper Root directory path.
But the stores are still created under the default location: $GEOSERVER_DATA_DIR/data.

I do not see any error message in the log files (most detailed verbosity level).

Has anyone experience with setting REST PathMapper Root directory path?


texts: Export coordinates as text from a shapefile layer in QGIS

texts:
How can I export X and Y coordinates from a shapefile layer as a text file in QGIS 3.18?
The height of every point is displayed in the attribute table, the x and y coordinates are not. They can only be accessed by using the object information tool and clicking on them in the map.


texts: QGIS 3.10 not loading layers from PostGIS

texts:
I am trying to load PostGIS layers in QGIS 3.10 but I don't load them, but in version 2.x I do load them normally. What will be happening, I attached an image of the error

No se pudieron obtener capas
Translated (es-ES to en-US): Could not obtain layers


texts: PostGIS: Extract intersected areas with their related column values

texts:
I have inserted a MultiPolygon table in PostGIS with the following columns:
id, e_id , p_id, com, geom

Assume that I have 4 polygons with overlapped areas with the following values:

id  com  geom
1   1_1  geom1
2   1_2  geom2
3   2_1  geom3
4   2_2  geom4

Is it possible to extract only the intersected areas and concatenate their related values ?
The result that I want to succeed is a table or shapefile like the following:

id  com      geom
1   1_1_2_1  geom1
2   2_1_2_2  geom2

If the polygon with id: 1 intersects with polygon with id: 3 and polygon with id: 3 intersects with polygon with id: 4

Until now I can extract only the geometries of intersected areas in a column with the following SQL code:

SELECT ST_INTERSECTION(a.geom, b.geom)
FROM mytable a, mytable b
WHERE a.ID < b.ID
AND ST_INTERSECTS(a.geom, b.geom);

Any suggestions?

UPDATE

with @JGH helpful answer, I managed to develop more my code and now I can concatenate the com values from two intersected areas.
The problem that remains the same is the id column that contains duplicate values. Is it possible to create a unique id column in the query ?

I tried SEQUENCE command in order to create a unique id column:

CREATE SEQUENCE serial START 1;
SELECT nextval('serial'),
a.com || '_' || b.com as com,
ST_INTERSECTION(a.geom, b.geom) as geom
FROM mytable a, mytable b
WHERE a.ID < b.ID
AND ST_INTERSECTS(a.geom, b.geom);

But it's not the solution that I wanted because, every time that I execute the code the id numbers increment.


texts:
This is related to another question I asked just recently:

Use OpenLayers WMSGetFeatureInfo on Mapserver Layer?

unicoletti indicated that he/she was able to use OpenLayer's WMS Layer declaration to generate and query MapServer WMS.  However, when I attempt the same, I get a error response:

mapserv(): Web application error. Query mode needs a point, imgxy and mapxy are not set.

This error makes sense in that when I use OL's MapServer Layer declaration, the url string includes these options with values filled in.  But using the WMS Layer declaration with standard WMS GetMap options, i.e. BBOX, WIDTH, HEIGHT, generates this error.  The image is not correct either which I describe next.

Another symptom is that the entire extent of the data defined in the mapfile is displayed rather than what's indicated in the BBOX when just doing a standard WMS GetMap request via a URL rather than using OpenLayers.

So for example, here's the first part of my mapfile:

MAP
  NAME "testgrib"
  IMAGETYPE      PNG
  EXTENT         -14000000 3000000 -7000000 7000000

And here's a GetMap request using standard WMS options:

http://localhost:8080/cgi-bin/mapserv?REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&FORMAT=image%2Fpng&MODE=map&MAP=wxmap%2Fglmpvis.map&LAYERS=mosaic&SRS=EPSG%3A3857&TRANSPARENT=TRUE&STYLES=&BBOX=-13807095.590531,3308267.5857033,-12241665.251469,4482340.3400001&WIDTH=640&HEIGHT=480

As you can see, the BBOX extent values are different from those in the mapfile.  But, MapServer returns an image with the extent values in the mapfile no matter what BBOX values are passed.  I've changed the mapfile extent to match the BBOX values and sure enough, an image of only that area is returned.

Hence this leads me to believe we've missed something during our MapServer installation, some setting or configuration.  I've checked to see that WMS was supported and received SUPPORTS=WMS_SERVER when making a GetCapabilities request as well as a -v on the MapServer executable.

I considered attaching the options to the query function that I can see are used to generate images using OL's MapServer Layer declaration to get around this issue.  But, as we eventually want our customers to access our data through the WMS, solving that would seem to be the main goal here.

I've searched the web and MS docs for a possible setting/configuration, but nothing really stood out as being critical.  I do have "wms_enable_request" "*" in the map file and, as indicated, I am receiving images from MapServer, albeit for the max extent of the data vs the viewing area.  And I'm receiving a response back to the query, so this doesn't look like a cross domain/proxy issue.  Also, this occurs with displaying an image for raster data as well as an image of point data.  The only thing that works is to use OL's MapServer declaration to generate the images properly for the viewing area.  But, that's not what we want for a long term solution with our customers accessing the data via WMS.

So, open to suggestions as to what we should either configure differently or what we might've missed to enable standard WMS in our MapServer.

Thanks!


texts: mapping over a million points using QGIS

texts: I'm trying to map a CSV with over a million lat/lon points in the US with no metadata. I've never tried to map this much data and want to know what steps I should take. I know it's too much data for QGIS so I tried splitting it into 100,000 line CSVs but then got a lot of files which I didn't know what to do with. I'm thinking it might be best to aggregate the points by census block but don't know the best way to do that using Javascript (I've never used Python) or command line tools? Any suggestions?


texts: moving user by gps and show interested points

texts: Is it possible in a Web application with google maps, to use a user his mobile gps and when move and is near on an interest point such as hotels, interested point i have stored in a Postgresql database, to show it in map.


texts: How to setup users and user groups on a distributed architecture?

texts:
I am learning ArcGIS distributed installation. I planned like this picture.

I read from installation documents about server users group.

agsAdmin
agsUsers

Which machine contains user groups and users? Web server, SOM or SOC PCs?
For example I installed Windows 2008 on  PC-1, 2, and 3 machines, all of these are default accounts. I installed PC-1 as SOM. PC-2 as SOC machine. Where my usergroups should be located or all pcs has these groups.


texts: How to import geometry data from csv file to postgres?

texts: I have geometry file in *.xls file, but when I save in *.csv it say "some feature in your workbook might be lost if you save it as CSV (Comma delimited)" however I keep that file anyway and when I import that file in postgres it say "ERROR: extra data after last expected column". Any idea please?


texts: Calculating Degree Days using iterate() and accumulate() functions

texts:
I am trying to calculate degree days using landsat's thermal band for a very specific geographic area and temporal range. I am having trouble coming up with a method that would allow me to get the result that I am looking for.

The degree day equation is quite simple: temp-base = degree units for that day. During a season (may-sep), farmers sum each day's total for a season total of degree days.

I want to be able to sum these degree days on a per-pixel basis, so what I imagine that would look like would be: pixel value - 285 = that pixel's DD. For every image that is in that range (every year has a different number of images), I want to sum up those values to give a total number of degree days for that pixel and then overlay that value onto the map.

This is the code that I have right now that finds the mean temp for a given range:

// Create a geometry representing an export region.
var features = ee.Geometry.Polygon(
        [[[-85.1269269325727, 44.839449231220215],
          [-85.2642560341352, 45.366750684688796],
          [-86.1926007606977, 44.92507820254033],
          [-86.4233136513227, 44.197108291897024]]]);
print("features",features);

// This function clips images to the ROI feature collection
var clipToCol = function(image){
  return image.clip(features);
};

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_8DAY_TOA')
                  .map(clipToCol);

var st = dataset.select('B10');

var imageVisParam = {"opacity":1,"bands":["B10"],"min":271.4,"max":284.6,"palette":["0f00ff","14b5ff","d3ff86","2aff5e","20c749","ffbd1e","ff931c","ff5018","ff0000"]};

Map.addLayer(
    st.filterDate("2014-04-01", "2014-06-01"), imageVisParam,
    '2014');
Map.setCenter(-85.75976, 44.85257, 11);

var imageRGB = st.mean().visualize(imageVisParam);

print(imageRGB);
Map.addLayer(
  imageRGB)
// Export the image, specifying scale and region.

Export.image.toDrive({
  image: imageRGB,
  description: 'LST_415_515',
  scale: 5,
  region: features,
  maxPixels: 76000000000,
});

I'm looking at this resource to try to better understand the iterate() and accumulate() function, but am having trouble fitting it into the context of this specific deployment.


texts: Automatically Publish Composite Locator Service with Rebuilt Address Locators

texts:
I'm using a composite address locator with ArcGIS Server 10.3.  It is based on 6 other address locators.

I'm trying to find a quick way to rebuild all of the address locators and publish the composite locator service to ArcGIS Server.

Is this possible through the command line or with Python?


texts: How to recalculate the the length of a geography linestring from nearby geography point in POSTGIS?

texts:
I have a geography type linestring
*route_geography geography(LineString,4326)*, its a route from point A to point B. Now I have a point C near that line. I want to find the nearest point on that line, lets call it D. And then calculate the distance between D and B.

Should I use *ST_Line_Locate_Point()* to find the point D?

And if I do find the point D should I then use *ST_Line_Substring()* to split up the AB linestring to make a new DB linestring?

These two functions are for geometry type not geography that my data is in, so what should I do?


texts: How to Buffer with intersections in QGIS

texts:
I am buffering a layer of facilities (points). I will then use the buffer to clip and intersect another layer, containing resources.

I want to answer the question, how much of the resource lies within X distance of the facility?

Buffer facility -> clip resource so only areas within buffer remain -> intersect facility buffer and resource to show amount of resource within distance of facility.

Now, I am using ftools to create my buffer. The problem is, when the facilities are already within X distance of each other it does not create a full circular buffer, the polygons created do not overlap.

How can I create a buffer in QGIS where the polygons of the buffer are allowed to overlap each other (possible in ArcGIS)?


texts: Problem Installing pgRouting 2.0; libcgal8 dependency

texts:
I've been trying to install PostGIS and pgRouting on an AWS Ubuntu server (first tried 13.10, then reverted to 12.04, both 64bit "small" instance) for several days now but no matter what I do I keep getting hung up by the dependency libcgal8.

Among many other things, I've tried using exactly the steps @dkastl recommended in this post:

sudo add-apt-repository ppa:ubuntugis/ppa
sudo add-apt-repository ppa:georepublic/pgrouting
sudo apt-get update
sudo apt-get install postgresql-9.1-pgrouting postgresql-9.1-postgis-2.0

After executing the last command the system complained that

Package 'postgresql-9.1-postgis-2.0 is unavailable, but  is referred to by another package

Checking the ubuntugis ppa: I determined that the correct name is postgresql-9.1-postgis. Substituting that package name got me to the roadblock I've encountered all along:

The following packages have unmet dependencies:
postgresql-9.1-pgrouting : Depends : libcgal8 but is not installable
E: Unable to correct problems, you have held broken packages.

I've tried everything to get around this issue including the following:

Separating the postgis and pgrouting install commands and switching the order they were executed
Tried using $ sudo apt-get install postgresql-9.1-pgrouting --no-install-recommends but this didn't work either (still complained about unmet libcgal8 dependency)
Tried downloading libcgal8 and its dependencies and installing them directly, which was a lost cause
Using $ sudo apt-get -f install which was suggested by the system at some point, and seemed to do something but it didn't fix the problem
So many other things that I either can't remember or were too desperate to mention

Hoping one of you can point to something really silly that this newb is doing or not doing.
Thanks in advance for your help.


texts:
I have a bunch of layers stored in PostGIS that I want to merge.  I have spent quite some time investigating various approaches but I ket tripping over getting the GIDs to update properly.

In the end the only way that I found that worked was to:

load the layers into QGIS from PostGIS and then save them as shape files in a temp dir.
remove the layers from QGIS and reload them from the shape files.
delete the GID attribute from the layers and save the shape files again.
then use Vector->data managment->merge shape files to create a merged shape file
finally use SPIT to get the merged file back into PostGIS

The purpose of this 'question' is two fold, firstly to document something at works and secondly to ask if there is a better way of doing it.

Things that I discovered on this journey:

you can not paste features into a PostGIS layer because these don't get the appropriate magic in the GID field to allocate a new number from the sequence. (this is arguably a bug in QGIS ??)
Splitting feature in a PostGIS layer works as it should and the new feature picks up a GID when the layer is saved.


texts: How to move the last vertex till the last angle is squared?

texts:
I have problem using Orthogonal Digitizing button from CAD Tools plugin menu.
It moves the first point of the polygon by holding ctrl key pressed and pressing right mouse button till the last angle is squared.
Is it possible to change setting to move the last vertex till the last angle is squared?


texts: How do you create vector tiles to represet landmass and oceans

texts:
I plan to self-host my own vector tiles and use Mapbox JS to render them in the browser for a design project where I create minimal abstract looking maps of the world. Hosting on Mapbox is not possible due to tile size limitations.

I've been able to get most of the data I need out of OSM (roads etc) just fine and process it using Tippecanoe. But for land and water, I'm having problems. Tippecanoe takes 24 hours to process the raw data and turns several mbs of shapefiles into several gbs.

I'd like to know if there are alternative approaches to the one I'm taking which will require less processing time.

Here is how I'm going about it.

For my use case I simply want to be able to render land as green and water as blue.

I do not need to distinguish between residential or industrial land. Nor do I need to know the difference between a river, lake, sea, or ocean. If it's land I want to color it one shade of green and if its water I want to color it one shade of blue.

I have downloaded Natural Earths 10m land and ocean files but they make Tippecanoe crawl. I think this happens because I need all my tiles to support a max zoom of 12 so I can utilize over zooming. I can't specify a max zoom of 6 for oceans only as if I zoom to 7 the ocean data will not show as the roads in the mbtiles have a max zoom of 12.

It seems rather inefficient to have to process 2gb of data for my use case. I'm wondering if I'm taking the right approach or overlooking a more simple solution?

One thing I considered was not to bother creating ocean data at all in the vector tiles. Instead I'll just color the "background" of the whole world map as blue (assuming this is possible in mapbox js). I'm not sure if this has shortcomings however.


texts: Filtering out OpenStreetMap data made by specific users using osmosis?

texts:
I am trying to find the code for a query in osmosis that will allow me to download only contributions made by certain OpenStreetMap users. I was thinking to use osmosis to filter out all the contributions made by some friends and I, in order to have some metrics on how many edits we made on OSM for a certain project.

Has anyone done this in osmosis?

I think the code should look something like this in osmosis, however this particular code is not working:

osmosis --rbf in liberia-latest.osm.pbf --nkv keyValueList="username.JohnDoe1" --wx ourcontributions.osm


texts: LiDAR dataset which is benchmarking for testing classification algorithms?

texts:
Is there a publicly available benchmarking database of classified (labeled) LiDAR data? For example, where each point is classified as ground, vegetation, building, road, water, etc?

The point clouds that I find listed in other questions are not complete classified (or complete validated) point clouds. I have yet to find one where each point contains a label that specifies the type of object from which it originated.

A set of this nature would be very useful as a way to quantitatively measure range finding algorithms. The scholarly articles that I have read typically rely on small private sets or upon qualitative measures.

I would be happy with any source of data, but I prefer airborne LiDAR. Also, I would like a broad range of data including cities with buildings, house, and trees and also pure forests, but again, I would take anything I could get.

In other words, I am looking for a LiDAR point cloud which all points were classified with more than 99.9% accuracy (probably manually classified point clouds?).


texts: Coding of data font non editable in shapefile in QGIS

texts:
What is the reason for a shape layer to have the field "Codificación de la fuente de datos" (coding of data font) selectable or not?

In the following image the red arrow shows the dropdown list in grey and changing it is not possible.

Sorry, my version of QGIS 2.18.11 is in Spanish.

I usually go to properties of the layer > General > Layer information, and there I can change the data font codification.

For more info: the shapefile is a real shapefile, it is possible to edit it, and it is saved in the same folder as other shapefiles that I can actually change this property from the drop-down list. Also my coworkers open it and can change the coding of the data font.


texts: Repeatedly losing data sources for raster layers in ArcGIS for Desktop?

texts:
I have several raster layers that keep losing their data connection whenever I add a new layer from running a tool. I repaired the data source, but this keeps happening every time I run a tool. Is there an option/preference that I need to change so this stops happening? I didn't move any files so the path names are the same.

UPDATE:
I restarted Arc and the document, which seemed to work for a little bit, but after a while it started happening again.


texts: Invalid Input Syntax for Type Double Precision in PostGIS?

texts:
I have this query that is trying to update my geom column that's built off of lat and lng coordinates

UPDATE main_table SET geom = ST_SetSRID(ST_MakePoint(lng, lat) ,4326);

However it returns that function st_makepoint(text, text) does not exist and gives me a hint to add explicit type casts, so I updated my query to

UPDATE main_table SET geom = ST_SetSRID(
ST_MakePoint(lng:: double precision, lat:: double precision) ,4326);

And now it's throwing an error invalid input syntax for type double precision:.

Just a note, right now my lats and longs are stored as text because I had to use a pg restore to load the table, so I'm sure that's where the problem is but don't know how to address it.


texts: How do we know the end of a street in OSM?

texts:
Is there a possibility to know the end of a street when I look in a osm-file? Is there a tag or a system which clearly identify the node which is the end/beginning? With the same street I mean here a street which have a connection and has the same name.

Some times there are streets with forks, which consist out of several ways - like 'Am Teich' in Greifswald (Germany).

But Often we have also several ways even for a street without forks, like Rue Billerey.

To illustrate my question, I made this drawing where we have three times the same points but three times different end-points of the street.


texts:
I am doing an analysis using census tracts from São Paulo, Brazil, with Qgis. The problem is that changes in their boundaries occur between the periods. Most of them divide in two or three. The literature usually says that we have to use a something called "minimum comparable areas", creating the smaller area possible and adding the attributes of the tracts that are inside this area to this new one. If the census tract remains unchanged between the periods, by example, this area corresponds to the minimum comparable area, and the attributes are his own.

Is there a way to do that procedure with two or more layers more automatically than one by one?


texts: In QGIS 2.18 I can view OSM using XY tiles no problem. In 3.2 however, copy/pasting the same string in the URL for a new connection doesn't work. The CRS for both is the same. In 3.2 it doesn't matter what I use, quickmapservices or xytiles, I cannot view tiles. Not Bing, Google, OSM, nothing. I must be missing something basic in the differences between QGIS versions but I can't think what it is.


texts: Minidump with Qgis Raster Calculator in python console

texts:
I want to calculate a new raster starting from my own raster. I use QGIS Raster Calculator and it works.
Then I want to include it in a python script to create it manually. My new raster is created but empty. I am searching where the problem comes from. My formula is correct (when I copy paste it in QGIS Raster Calculator it works fine).
So I tried to execute it directly from the console, I got this error and qgis closed :

My code is the following :

l=iface.activeLayer()
a=l.extent()
b=l.width()
c=l.height()
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
rast = QgsRasterCalculatorEntry()
entries=[rast]
calc = QgsRasterCalculator( '("lu_6@1"=210)*200',"D:/ju/Histo_Europe_folder_test/test_map.tif",'GTiff',a,b,c,entries )

My raster layer is called lu_6 and I am using QGIS 2.18.9.

Any idea why I get this error ?


texts: How to add GDAL Programs folder to ~/.bash_profile

texts:
I have recently installed on Mac OS X Lion 10.7.4 the GDAL Framework 1.9 with Python bindings from [http://www.kyngchaos.com/software:frameworks][1].

I need to add the following Programs folder to my ~/.bash_profile as described in the ReadMe file:

export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH

This is what's currently on my ~/.bash_profile:

# Setting PATH for Python 2.7
# The orginal version is saved in .bash_profile.pysave
PATH="/Library/Frameworks/Python.framework/Versions/2.7/bin:${PATH}"
export PATH

What is the right syntax for adding the new path?

I am new with everything (Mac, Python and GDAL) so please excuse me if I wasn't clear.

Many Thanks

Andrea


texts: Combining continuous raster and point in ArcGIS for Desktop?

texts:
I have a layer that's a topographic map of a mountain range, represented by a continuous raster. I have another layer that is just points, representing where the summits are.

How can I use the raster layer to attach an elevation to each of these summit points based on location?

I'm using ArcMap 10.1.

I've tried using Extract value to points but it will not work, saying "Failed to add RASTERVALU to table".


texts: Get majority value of raster data inside polygon using PostGIS

texts:
We have imported a tiff with landcover data into PostGIS as a raster table.
We also have a polygon table with provinces.
I'm trying to get the majority of the landcover in each province.
But can't get it to work. I already looked at Return Majority from PostGIS function ST_SummaryStatsAgg() and we are using ST_SummaryStatsAgg() already with other raster data. But we are stuck with the Majority calculation.
The problem we have with the mentioned link is that we don't have the class attribute nor do we know how to create it from the raster data.


texts: Tracing path of lines in network using QGIS

texts:
I am using QGIS 3.

I want to trace the path of lines in network. I found that I can use Snapping Toolbar to trace the lines but in manual and I want to do it in automatic, I give start node and it automatically traces the path.

I have a map and I want to trace the path for each line from start node till the end, then extract each line id. But I want to it in automatic method, not in manual method as snapping tool do.


texts: Rotating around arbitrary point, not center when using Google Earth Pro to Edit Image Overlay?

texts:
In Google Earth Pro "Edit Image Overlay", we can rotate the overlay
around its center (via grabbing the green diamond,) but how might one
rotate it around an arbitrary point not its center?

(Yes, doing so would just be a shortcut for a rotation around its
center combined with a translation (grabbing the "+" at the center and
moving it.))


texts: Exporting attribute table - some records missing

texts:
I'm using ArcGIS 10.1 (SP1 for Desktop) and when I'm exporting my attribute table to a dBase  format, it will export some but not all the rows. There is no error message indicating anything has gone wrong. I only notice it when I'm opening the table in Excel. I've been experiencing this issue with several different tables (each of varying sizes from 10 records to 300000) and this problem doesn't happen every time. The example I'm going to show is with one of my larger tables.

Here are my steps for exporting an attribute table (with my selection already completed):
Table Options > Export > Selected Records > Choose Location > Name > dBase Table > Export

This is my attribute table. I need to export 174230 rows.

My exported table gave me 51125 rows. This is what my exported table looks like:

After my several tries, I was finally able to export the full table. I had to put a definition query on. Sometimes, I get around the problem by doing the standard copy and paste. This depends on the size of the table I'm working with.

I know there are workarounds, but this problem is quite persistent. Has anyone experienced this problem? Is there something wrong with my settings?


texts: Extent from Centroids, Projection, and Resolution?

texts:
I have a file that defines a raster, although all the metadata was lost.  I know the lat/long coordinates of the grid centroids, and that the grid is in a Mercator projection with a 1/12deg equatorial resolution.  Is there any way to define the extent of the raster by just knowing the centroids and resolution?  Currently I'm doing:

xmin = min(longitude)-1/24
xmax = max(longitude)+1/24
ymin = min(latitude)-1/24
ymax = max(latitude)+1/24


texts: ST_TRANSFORM SDE Oracle does not convert the geometry

texts:
I try to convert a geometry from 31370 to 3857. But it does not work. The coordinates stay in 31370.

DB: Oracle 11

ArcSDE: 10.4.1

The code:

select
sde.st_srid(shape)
,sde.st_x(shape) as X_lambert
,sde.st_y(shape) as Y_lambert
,sde.st_x(sde.st_transform(shape, 31370, 3857 ) ) as X_pmercator
,sde.st_y(sde.st_transform(shape, 31370, 3857 ) ) as Y_pmercator
from usersde.cst
;

The datas:


texts: Using Update Cursor in a double for loop in arcpy

texts:
I have some problems using a double for loop and getting results written into a geodatabase.

I have got a shapefile with polygons (param_polygone) and I want to add data to it, which is stored in a table (table_datei). In the polygon shapefile I have got features with an unique ID, called "FID". In the table I have got fields, which store lenghts ("LENGTH") belonging to a polygon (src_FID). So often I there are several lengths stored for one polygon. For example there is one FID in the polygon shapefile, which is 0 and several entries or rows in the table, where the src_FID is 0 and where unique lengths are stored in the same row. I want to sum up all lengths of one src_FID and store it in the field "Nbr_Laenge", which belongs to the FID that is the same value as the src_FID.

Somehow only the first value gets sumed up in "Nbr_Laenge". All the others have the value 0 in the end of the execution.

Here is what I got so far:

table_cursor = arcpy.da.UpdateCursor(table_datei, ["src_FID", "LENGTH"])
geb_cursor = arcpy.da.UpdateCursor(param_polygone, ["FID", "Nbr_Laenge"])
    for geb_row in geb_cursor:
      for table_row in table_cursor:
        if table_row[0] == geb_row[0]:
          geb_row[1]=geb_row[1]+table_row[1]
          geb_cursor.updateRow(geb_row)

Something with the for-loops isn't working properly. I am new to python and working with ArcGIS.


texts: Leaflet Routing Machine: how to dinamically change router settings?

texts:
I'm using Leaflet Routing Machine in a litte Leaflet + Javascript map.

I'm using Graphhopper as router in this way ....

var control = L.Routing.control({
  waypoints: [
    L.latLng(57.74, 11.94),
    L.latLng(57.6792, 11.949)
  ],
  routeWhileDragging: true,
  geocoder: L.Control.Geocoder.nominatim(),
  router: L.Routing.graphHopper('a585904f-5193-4605-bc3c-870c4f472177' , {
    urlParameters: {
        vehicle: 'car'
    }
  })
}).addTo(map);

..... and all works fine

Now in my map I'd like to put give the opportunity to select, using a couple of buttons, the router "vehicle type", so if the user select "foot", I'd like to change the  route vehicle parameter setting to 'foot' value.

How may I set this value? Something like

 control.Router.?????? = 'foot'


texts:
I have line like Road and points are like housenumbers!

I have fields like side and ID_T

Highlighted blue apts having same ID_T because hooking one edge
and side, if it hook right get R or else get L

I need to select which apts are hooking only one side


texts: QTest interactions with the QGIS Map canvas

texts:
Is it possible to simulate right or left mouse button clicks on a QGIS map canvas.  Currently writing tests for a plugin and need to be able to create a line geometry thorough the use of theses buttons.

Have been trying...

QTest.mouseClick(canvas, Qt.RightButton, pos=QPoint(1200,500), delay=1000)


texts:
I tried to add several Google Maps Engine layers using the following code:

gmeLayer = new google.maps.visualization.DynamicMapsEngineLayer({
    accessToken: token,
    layerId: layerId
});
gmeLayer.setMap(map);

I would like to display the ploygon layer on the bottom and point layer on the top.

How can I set the order for each layer?


texts:
I want to get GeoJSON data for the UK. The regions which I require are in the map displayed in the below link:
Former Government Office regions

I have been searching the net but all in vain. Can anyone help?


texts:
Every time I execute this, I get a runtime error pointing to my "for row in cursor," stating that A column was specified that does not exist. Any ideas how to fix this?

import arcpy
import os

gdb = r'C:\Users\Stephanie\Documents\SPRING20\PGIS\project1\usa_cities.gdb'
feature_class = "largest_cities"

fc = os.path.join(gdb, feature_class)
print(fc)

fields = ['NAME', 'POP2017']

with arcpy.da.SearchCursor(fc, fields) as cursor:
    for row in cursor:
        if row[2] != None:
            print("City: " + str(row[2]) + " Population: " + str(row[9]))


texts: Warning 000405 when trying to calculate field

texts:
I am trying to create a model but am getting hung up on my first step. I need to calculate an ID field from an FID field so that I can join it to another table based on the ID field. So I have a table:

FID  Shape     ID
0    Polyline  0
1    Polyline  0
2    Polyline  0

Through the calculation, I was hoping to get something like this:

FID  Shape     ID
0    Polyline  1
1    Polyline  2
2    Polyline  3

Easy expression: [FID]+1. When I run the calculation though, I get the error 000405: No records within table. Since there won't be any data in the table until I join it to the other table, I'm not sure what to do here. Even trying to create a new field and running the same calculation returns the same warning.


texts: Workflow process to style vector data using GeoServer

texts:
I'm starting to get a better understanding of map-servers in general, but I am still unclear on how best to use GeoServer. At the moment, I use web-hosting to store the web-browser related parts of my site, and php code to connect to a postgresql/postgis database. All of my data is stored in a desktop server, which I chose to use as it is a very low cost way of storing a large amount of data.

1] web-browser <--> database

This works fine also, as I do not need to style vector data; I am only returning numerical results with geojson vector overlays. However, I would like to incorporate GeoServer into this workflow so that I can apply better cartographic styles to data in the future. I understand the work flow should be similar to the following:

2] web-browser <-->  GeoServer <--> database

My difficulty in understanding this piece, is how to connect GeoServer into the process. Can I run GeoServer on my desktop and connect to it to the browser via php, similar to the database approach? Or does GeoServer need to be installed on the web-host and managed there?


texts:
I intend to use QGIS 3.6 as a tool for archaeological field works, live tracking.

I have created a form and everything works great except for creating polygons. I have installed QGIS on a Panasonic FZ-G1, Windows 10, and use an external GPS. Points and lines work well with GPS.

When I try to create polygons with GPS, everything works well until I finish the polygon, then a break point ends at the equator, which is completely wrong. If you then validate the polygon, the result is that it is not closed. I have tried several older versions of QGIS, I have tried different GPS receivers and tried on different computers and with different projections. The error is the same.


texts:
I want to get some statistics about my geonetwork catalog. Things I would like to know are

total dataset count
dataset count for different format types (e.g. CSV, WMS, etc.)
dataset count that is valid for a specific geographic region

The best thing I could think of so far, is to query the catalog with the GetRecords operation, optionally filtered to restrict feature count.

Then it would be neccessary to parse the response and count the occurrence of e.g. individual identifiers.

This would be possible but doesn't sound too nice, as I assume there should be some better solution already available for this, because it sounds like a must have for a geodata catalogue.

Does anyone know of a better solution to this, by means of an OGC request to csw or a geonetwork specific solution?

UPDATE:

Just found, there is a built-in solution to this in the catalog UI as described here. Thats actually what I need, but I would like to query it from the API to show it in tabular format.


texts: How to change the color and width of the shapefiles via arcpy

texts:
I have many Shapefiles in the same folder. Each shapefile contains polylines. I want to change the color and width of the lines. I'm quite new to ArcGIS and hope someone can help me.

  ApplySymbologyFromLayer_management (in_layer, in_symbology_layer)

Is this what I have to use? How to generate an "in_symbology_layer"?

The following code is what I have so far and I will edit it if I know any information.

import glob, os
import arcpy

the_path = r"C:/project/"
os.chdir(the_path)

for file in glob.glob("*.shp"):

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

    df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

    newlayer = arcpy.mapping.Layer(the_path+file)

    arcpy.mapping.AddLayer(df, newlayer,"TOP")

    arcpy.RefreshActiveView()
    arcpy.RefreshTOC()

    del mxd, df, newlayer


texts: Custom c# GP Tool ArcGIS Server 10.1 - build for 'Any CPU', x86, or 64bit? Does it matter?

texts:
I think I understand most explanations of these targets, but I get confused in the context of ArcGIS customizations like this because what we are compiling has to work in conjunction with ArcGIS, not just by itself.

We are having trouble getting our tool to register on ArcGIS Server 10.1 with esriregasm, though it works on Desktop. Build platform is one variable, so I thought I'd ask which I should build for when the intent is for the tool to be installed on ArcGIS Server 10.1.

Thank you!


texts: Using Postgis,The shell prompt Invalid endian flag value encountered,when I want to do MBR to my WKB_geometry

texts:
I want to calculate the MBR for the column of table.

I have already install Postgis and postgresql.

the data is cvs

and I load the data to my postgresql-table by [copy]

and my table scheme is

and When I running

  SELECT ST_Envelope (wkb_geometry::geometry) from test1;

OR
running

  SELECT PostGIS_ADDBBox(wek_geometry) from test1;

the shell prompt me

  ERROR:Invalid endian flag value encountered.

I just want to create the Minimum bounding box（MBR） to my table(by postgis,and using web_geometry column),but I don't know what should I using the function and what the correct syntax


texts: How to configure wampserver for spatialite

texts:
I am in windows 7. I installed WAMPServer. Now, I can not load spatialite libraries. It is showing warning..

  Warning: SQLite3::loadExtension() [sqlite3.loadextension]: Not
  supported in multithreaded Web servers

Here is my total configuration procedure, what I have done...

From this link I got spatialite lib. I copied "libspatialite-1.dll" and paste it to "D:\wamp\bin\php\php5.3.8\ext" which contains php extention dlls.

http://www.gaia-gis.it/spatialite-2.3.1/libspatialite-win-x86-2.3.1.zip

Then I edited php.ini file. I changed the following configuration.

.
sqlite3.extension_dir = C:\libspatialite-win-x86-2.3.1\bin
.
extension=libspatialite-1.dll
.
enable_dl = On

And lastly, I copy and paste all the libraries from my downloaded libspatialite-win-x86-2.3.1 to my project folder(libspatialite.a, libspatialite.dll.a, libspatialite.la) in my php code i write the script as follows...

<?php

$db = new SQLite3('sixcommunes.sqlite');

$db->loadExtension('libspatialite.a');

$rs = $db->query('SELECT spatialite_version()');
while($row = $rs->fetchArray()){
    print "
SQLite version: $row[0]
";
}

?>

I do not know what I have done wrong or how to solve this problem?


texts:
I'm trying to do intersect with 2 shapes, but it's not working because it can't read my shapefiles... I don't know why, because when i ran it with my other script, it works perfectly!

Traceback (most recent call last):
  File "Z:\04 - Temp\Charles\Teste_python\new toolbox.tbx#Script_newtoolbox.py", line 106, in
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 335, in Intersect
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 332, in Intersect
    retval = convertArcObjectToPythonObject(gp.Intersect_analysis(*gp_fixargs((in_features, out_feature_class, join_attributes, cluster_tolerance, output_type), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing_base.py", line 506, in
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
 ERROR 000732: Input Features: Dataset 'Z:\04 - Temp\Charles\NDVI_tempo\NAL\Zonalpoly.shp' #;MA02ASEZ3S.shp # does not exist or is not supported
Failed to execute (Intersect).
Failed to execute (Script).

NDVI_tiff = Divide (Rastdiv1, Rastdiv2)

            Per_tiff = Times (NDVI_tiff, 100)

            fc_rast = 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\fc_rast.tif'

            Zonal_tif = ZonalStatistics(fc_rast, "ID_TALHAO", Per_tiff, "MEAN", "DATA")

            Zonalint = Int (Zonal_tif)

            arcpy.RasterToPolygon_conversion (Zonalint, 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\Zonalpoly', "SIMPLIFY", "VALUE")

            Zonalpoly = 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\Zonalpoly.shp'

            inFeatures = [Zonalpoly, fc]

            arcpy.Intersect_analysis(inFeatures, 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\intershp.shp', "","", "")

This is the script, where "fc" is each shapefile in my workspace... It reads the name of the file (MA02ASEZ3S.shp), but the error 000732 seems to warn that something is wrong...

And this is how the fc variable is declared

down = 'Z:\04 - Temp\Charles\NDVI_tempo\NAL_PL'
image_server = 'Y:\E - IMAGEM SERVER\IMAGEM_SERVER.gdb\F02_RAPIDEYE'

arcpy.env.workspace = 'Z:\04 - Temp\Charles\NDVI_tempo\NAL_PL'
fcs = arcpy.ListFeatureClasses()

for fc in fcs:

    arcpy.DownloadRasters_management (image_server, down, "", fc, "CLIPPING", "CONVERT_AS_REQUIRED", "TIFF", "NONE", "", "NO_MAINTAIN_FOLDER")

    arcpy.PolygonToRaster_conversion (fc, "ID_TALHAO", 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\fc_rast.tif')


texts: Why does thumbnail on one MXD not update?

texts:
I want to make modification on multi-mxd properties, I try Arcpy script with a loop like " for mxd in mxdList:" , but I always have same results, Thumbnails are not created...

import arcpy,os
dossierRech = "C:\TEMP"
arcpy.env.workspace = dossierRech
mxdList = arcpy.ListFiles("*.mxd")
for mxd in mxdList:
    filePath = os.path.join(dossierRech, file)
    mxd = arcpy.mapping.MapDocument(filePath)
    mxd.deleteThumbnail()
    mxd.makeThumbnail()
    mxd.save()

I try on only one MXD , it doesn't created it too.

Is there some issue with "mxd.makeThumbnail()" ?


texts: ArcGIS 10.1 ModelBuilder: Data displayed, Multiple ring buffer causes Error 999999

texts:
I am trying to make a flow of processes in ModelBuilder, connecting multiple tools:

import polygon dissolved it,
intersected with another files, and
selected desired polygons by attributes("ZONA" = 'A' AND
"Shape_Area" >= 30).

UNTIL HERE WORKS GREAT
create buffers by Multiple Ring Buffer
            - DOES NOT WORK !!!!

I can successfully display my Selected polygons (OUTPUT).

However, when I run "Buffer" or "Multiple Ring Buffer" tools in ModelBuilder on OUTPUT , I've got back: 999999 Error executing function.

What can be wrong with my OUTPUT data, if I can display them, but can't run the buffer analysis? And more important, how can I fix that error? Maybe save my polygons by another tool (something like Collect values ..) that would prevent the problems?

EDIT:

All my data have defined coordinate system, are single part polygons (results from Multipart to single part Tool), cover area max of 6000 hectares (60 km2).
I am working on ArcGIS 10.1, Windows 7 Ultimate 32 bit.


texts: Extracting shapefile with suburbs of a city from OSM

texts:
I never downloaded shapefiles from OpenStreetMap before and now I am trying to download a shapefile of the city of Turin (Italy) with the suburbs polygons encoded in it. I know for sure that the information about suburbs is there in OSM (check here), but I could not manage in any way, using any tool, to download something different from the circoscrizioni areas, which are bigger. If I download the Turin shapefiles from, for example, websites like bbbike I get the usual shapefiles (buldings/landuse/natural/points etc.) but the suburbs are nowhere to be found.

Any suggestions?


texts: Numerical fields decimals limited by default

texts:
I'm using ArcMap 10.5.1 and working from a geodatabase.

I just want to limit the number of decimal numbers shown to 1. Currently 6 decimals are shown, and I know I can reduce them manually by opening each field properties window and setting the value to the desired decimals.

The problem is that I'm working on a ModelBuilder to get the desired result for my project, and I have to run it every once in a while, therefore losing the manually edited fields.

Also, the value I'm interested in this project is the polygons area, which is an automatically created field and I can't setup it's properties from the beginning.

Is there any way to set up how many decimals will be shown by default in the field everytime a new one is created?


texts: Combining single cloud-masked images using Google Earth Engine

texts:
I am trying to calculate the mean pixel values of three Sentinel-2 images, all of which I have cloud masked separately. They all have small areas of cloud but in different places. I assumed that when I calculated the mean, the cloud masked pixels having no-data values would simply be excluded from the analysis and in cloudy patches I would get a mean of the other two images instead. However this is not the case and my final image is masked wherever there was cloud in any of the three images.
My goal is therefore to combine all three of these images, masked for cloud, to produce a cloud-free image of average pixel values.
I considered that calculating the mean directly (instead of adding separately and dividing by 3) might solve this, however since my images are separate assets and not a filtered image collection I do not know how to do this.
function maskS2clouds(collection) {
  var qa = collection.select('QA60');
  var cloudBitMask = 1 <<10;
  var cirrusBitMask = 1 <<11;

  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
  .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return collection.updateMask(mask).divide(10000);
}

//image collection is filtered
var dataset = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .filterDate('2017-06-19', '2017-06-20')
  .map(maskS2clouds);

var dataset2 = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .filterDate('2017-08-23', '2017-08-24')
  .map(maskS2clouds);

var dataset3 = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(geometry)
  .filterDate('2017-07-19', '2017-07-20')
  .map(maskS2clouds);

//image with mask applied
var image1 = dataset.map(function(image) { return image.clip(geometry); });
var image2 = dataset2.map(function(image) { return image.clip(geometry); });
var image3 = dataset3.map(function(image) { return image.clip(geometry); });

var june = image1.mosaic();
var aug = image2.mosaic();
var july = image3.mosaic();
var junaug = june.add(aug);
var augjuly = junaug.add(july);
var avg = augjuly.divide(3);


texts:
I have made a map using Mapbox Studio. The map is made up of 2 map sources,

MapBox Streets
Custom Style - a PostGIS layer

I am having trouble working out how to set the draw order. I have tried the logic of setting it with the CartoCSS and also tried dragging layer "Data" down the list (image below) to affect the draw order and have had no luck because I can't drag the layers.

Any points from the Geo Ecosystem?


texts: Labeling large number of streets in ArcGIS

texts:
Anyone know of a decent way to label a large street network for large scale printing that doesn't involve manual editing/shifting? Even with advance placement properties there seems to be a number of overlaps and just general poor placement of the labels. This is not a straight forward network as there is a ton of variation between road lengths and shape.

Any plugins exist that might help or other programs that do a better job at labeling?


texts:
I have a large number of polygons, and I'd like to calculate an approximation of the distance from each one to each other one.

However, calculating the distance between two polygons is slow.  Meanwhile, approximating each polygon as its centroid and then calculating distance is very fast, but not so accurate.  What approximation of a polygon can I use so that I can obtain both quick and (reasonably) accurate distances?


texts: Import own modules in custom forms in PyQGIS

texts:
I want to create my own custom forms using Python for QGIS, and since some of the forms will share some definitions and objects I would like to have my own classes for all of them. For that, when I try to import my classes in the python logic form file I get this error with the following tries:

With import comDef I get ImportError: No module named 'comDef'
With from .comDef import *  I get SystemError: Parent module '' not loaded, cannot perform relative import
With from comDef import *  I get ImportError: No module named 'comDef'

These custom forms are located in a folder where the custom python logics and custom classes are.

Can you provide any hint about this issue?


texts: ArcGIS Modelbuilder Warning empty output generated to Halt process?

texts:
Is there a simple way of making ArcGIS stop if the following message comes up after "copy rows" have finished (see picture);

  Executing (Copy Rows (4)): CopyRows QueryTableB_cx
  C:\Users\fmromss\Documents\ArcGIS\Default.gdb\QueryTable_CopyRows2_cx
  Start Time: Mon Apr 25 12:40:25 2016 WARNING 000117: Warning empty output generated. Succeeded at Mon Apr 25 12:40:25 2016 (Elapsed
  Time: 0,38 seconds)

It would be quite neat if It did not export it to excel if theres not any actual data.

Bear in mind this is just a small part of the model hence I dont want the entire model to stop, just this process.


texts: Is there any way to use open source tools to make sure that lines on a map don't overlap?

texts: I am a generalist developer getting the hang of basic GIS. I made this GeoJSON map by clipping and converting different shapefiles. The problem is that the precincts of New Orleans extend into the Mississippi river. I want them to stop at the river. Is there any easy tool or method to get the shapes from not overlapping? The borders of the precincts should end at the river. The precincts should not extend into the river. I am using all open source tools.


texts: Finding Closest Point to Polyline Vertex and Apply Elevation

texts:
I am trying to code a tool that will evaluate all vertices along a polyline and then find the closest point from a survey shapefile and apply that elevation to that vertex, I am wanting to iterate through a bag of vertices and then at each vertex iterate through the points shapefile and use the distance between the points and elevation in a dictionary and then find the smallest distance at each vertex and grab that elevation here is the code I have so far.

Here is a map package with a typical map https://drive.google.com/file/d/12jrnIn6BAaduLx10_6sI2OJhWB7-IBr2/view?usp=sharing

import os
import arcpy
from math import sqrt
##Declare map and layers, split layer from all layers and then parse .shp    files
##declare map
mxd=arcpy.mapping.MapDocument("CURRENT")
##dataframe
df=arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
 ##pipelayer
listLayers=arcpy.mapping.ListLayers(mxd,"PL_PIPELINE_LN",df)
staging="c:\\staging"
fileName="Report.txt"
path=os.path.join(staging,fileName)
f=open(path,"w+")
##point features to get elevation from
where=""""CODE" = 'BEND' OR "CODE" = 'FLN' OR "CODE" = 'NIP' OR "CODE" = 'RED' OR "CODE" = 'TAP' OR "CODE" = 'TEE' OR "CODE" = 'VLV' OR "CODE" = 'WELD'"""

#lists and dicts for values
pointDict={}
pointList=[]
pointPair=[]
featureList=["FLN","NIP","RED","TAP","TEE","VLV","WELD"]
featurePoints={}
lines=[]

#spatial reference to measure stateplane values
spatialRef=sr_state=arcpy.SpatialReference(6578)

#get centerline name so we can select those features
cenName=""
with arcpy.da.SearchCursor("PL_PIPELINE_LN","CENTERLINE_NAME") as sc:
    for row in sc:
        cenName = row[0]

##get count of features for iteration  , haven't used this yet but may to control iteration
arcpy.SelectLayerByAttribute_management("PL_PIPELINE_LN","NEW_SELECTION","""CENTERLINE_NAME = '{}'""".format(cenName))
result=arcpy.GetCount_management("PL_PIPELINE_LN")
count=int(result.getOutput(0))
print(count)

##clear selection and begin process
arcpy.SelectLayerByAttribute_management("PL_PIPELINE_LN","CLEAR_SELECTION")

##select line features with cen name
with arcpy.da.SearchCursor("PL_PIPELINE_LN",["SHAPE@","OBJECTID"],spatial_reference=spatialRef) as sc:
    for row in sc:
        ##access geometry and get into vertex
        geometry=row[0]
        objId=row[1]
        print (objId)
        pts=geometry.getPart(0)
        for point in pts:
            ##for vertex in vertices
            #print (pointList)
            linePoint=arcpy.PointGeometry(point)
            #featurePoints={}
            ##clear pointDict for each vertex
            pointDict={}
            ##iterate through the survey shape file store values to point dictionary
            with arcpy.da.SearchCursor("1708031_20180611",["SHAPE@","CODE","ELEV","FID"],where,spatial_reference=spatialRef) as sc2:
                for row2 in sc2:
                    pointGeometry= row2[0]
                    code=row2[1]
                    elev=row2[2]
                    oid=row2[3]
                    shpPoint=pointGeometry.getPart(0)
                    ##get distance to vertex
                    dist=sqrt((point.X-shpPoint.X)**2+(point.Y-shpPoint.Y)**2)
                    ##print(dist)
                    ##add values to point dictionary
                    pointDict[oid]=dist,elev,code

            ##find smallest distance in dictioanry
            minDist=(min(pointDict.items(), key=lambda x:x[1]))
            print(minDist[1][0])

        ##minDist=(min(featurePoints.items(),key=lambda x:x[1]))

I get good returns but I don't always get a minimum I suspect my dictionaries aren't complete.


texts: Editing an ArcGIS Online Hosted Feature Service with ArcGIS Desktop Basic

texts:
I'm working to develop a workflow to edit ArcGIS Online hosted feature services on desktop. The Create Local Copy for Editing and Synchronize Local Edits with Server functionality works as expected for my ArcGIS Desktop Advanced license, but the client only has ArcGIS Desktop Basic.

What's the recommended workflow in this situation? Our current approach is to Create Local Copy for Editing, then overwrite the existing feature service on ArcGIS Online. This seems a little heavy handed given how easy it is to edit this data on the web. Are there other options for Desktop editing with a Basic license?


texts: ArcGIS Online Edit/Smart Edit Widgets Not Working

texts:
I've configured a web app in ArcGIS Online and I'm unable to get the editor and smart editor widgets working.  Both widgets state the web map has no editable layers.  I've enabled editing for the layer I'm working with and I've confirmed I'm working with the correct web map.  I think the issue may be a bug in the software because I'm able to edit the layer in map viewer and I can edit the layer using the predefined edit template.  Could the issue be sharing?  I'm currently only sharing with my organization.  I guess it would help if somebody could confirm if the steps below are missing anything.  Again, it's odd because the web map does have editable layers.

Step 1: Upload layer to AGO

Step 2: Configure sharing with my organization

Step 3: Enable Editing

Step 4: Create Web Map, Add layer, (editing works at this point)

Step 5: Share web map

Step 6: Create app using the web app builder

Step 7: Add the editor widget (no editable layers at this point)


texts: Adding a grid and selecting random cells from spatial polygon in r

texts:
I have a fairly similar question to this: Similar Question
I am working with a map of Denmark or more a subset of a Denmark Map.

This file comes from the Global administrative areas site GADM site

require (raster)
require (rgeos)
require (rgdal)

DK1 <- getData('GADM', country='DNK', level=0)

#create an extent for clipping in raster package
ext=  extent (8.41, 9.04, 56.92, 57.16)
clipe <- as(ext, "SpatialPolygons")
proj4string(clipe) <- CRS(proj4string(DK1))
cropDK <- SpatialPolygonsDataFrame(clipe, data.frame(x = 1), match.ID = FALSE)

#clip the shapefile using rgeos package
DK_reg1<- gIntersection(DK1, cropDK)
plot (DK_reg1)

I would like to create a grid of 1000 equally sized cells, plot it and then select and highlight 100 random cells. I seem to be having some trouble because of the projection and undefined units. Does anyone have any tips on how to fix these issues?

#Grid the clipped shapefile to 1000 equal size units.
proj4string (DK_reg1)

bb <- bbox(DK_reg1)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration); 1 m = 3.28084 ft
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd

sp_grd <- SpatialGridDataFrame(grd,data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(DK_reg1)))

Finnally is it possible to add this gridded shape file to a ggmap raster showing the main roads?

library (ggmap)
library(mapproj)
map <- get_map(location=c(lon=8.65, lat=56.955665), zoom =10)
ggmap(map, extent = "normal")
  plot (DK_reg1, add=T)


texts: Selecting column like Primary Key in View from PostGIS using PyQGIS?

texts:
Postgresql doesn't support Primary Key in views.
When I load a postGIS view in QGIS "manually", I can choose the column that will be the Primary Key, but I would like to know if I can select a column as Primary Key when I load it using pyQGIS script.


texts:
I would like to store an r-tree or quadtree index in a sqlite database using python.

Does anyone know if there is any code to do this?

A bit of background:

The dataset cannot be moved to another format.
the dataset in question is a shapefile
Though there is spatial indexing for shapefiles, it's not meeting the needs
I want to avoid using spatialite extension (it's messy to install on windows)

I can create a quadtree/rtree in-memory, but since the dataset is very large, the index takes a longtime to build. In addition, it consumes a large about of memory and effects other analysis down stream.  If I could store the index in sqlite it would help get past that issue I think.


texts: How can I use polygon grid to export raster

texts: I have a polygon grid which consists of 100 square polygons. For each square I need to export the raster overlay. I have ArcGIS 10.0, FME 2014, Saga-GIS. Is there a suitable tool for this or do I need to script something with Python or Model builder?