texts: Creating new lines from existing lines based on points lying along the lines

texts:
I am pretty new to PostGIS and I have a problem I just can't figure how to start. I have a table of lines (id, name, the_geom) and a table of points (id, name, the_geom). Each point is located on a line and a line might have one or more points etc. The points can be located anywhere along the line, not just at vertices.

How can I create new lines from each line using points along the line given just the lines and points tables. For example, line ABCDEFG is a single line with its unique name and unique id in the lines table. How can I create new lines AB, BC, CD, DE, EF, FG such that they all have the same name as line ABCDEFG but different id, and load their new geometry into a new table of lines using PostGIS? Note that A and G are the start and end of line ABCDEFG and not actual points in the points table.

Example 1:

This is not an answer, but an attempt to explain the question better, following Chris' comment.

@Chris, the links you gave are nice but quite different from what I want to achieve. I am not trying to get the individual segments of lines based on the vertices of the lines themselves. What I am trying to achieve is to make sublines from lines based on where points from a different table intersect the lines. I think my diagram is a bit confusing, I have uploaded something better.


texts: Min and max values not correctly calculated for DEM using QGIS

texts:
I'm working with DEM in QGIS and need to know the min and max values of my data.
In the Properties window, under the Metadata tab, the statistics section shows such values. However, I usually find some pixels with values higher than the max and lower than the min.

Is there another way to obtain these values?

I've tried this with several DEMs.


texts: Detecting if process in model has successfully run in ModelBuilder?

texts:
I have a model designded in ModelBuilder to process raster datasets. It embraces a whole bunch of processes such as raster calcualtion, BlockStatistics and Kriging. I am wondering how I can put in success validation for each of the process in the model in ModelBuilder to detect if the process has successfully run.

Is it possible without converting the model to python?

The Model will be eventually published as Geoprocessing Service by ArcGIS Server 10.1 so that the web browse calls the service and then gets a Result as the output in JSON if the Model rans successfully or on which process it failed. The success validation component should be part of the Model and generate result as parameter of the model.


texts:
I have this problem, I need to paint each polygon with a distinct color based on an atribute value using .NET, for example:

Polygon    value
1        10
2        10
3        5
4        5
5        1

Polygon 1 and 2 will be painted blue, polygon 3 and 4 will be green, and polygon 5 will be painted red

Does anyone knows how to do this? Thanks for the help.


texts:
When I load OpenStreetMap data to PostGIS using osm2pgsql the resulting tables look like this:

+----------+--------------------+--------+---------+
| Schema   | Name               | Type   | Owner   |
|----------+--------------------+--------+---------|
| public   | planet_osm_line    | table  | user    |
| public   | planet_osm_nodes   | table  | user    |
| public   | planet_osm_point   | table  | user    |
| public   | planet_osm_polygon | table  | user    |
| public   | planet_osm_rels    | table  | user    |
| public   | planet_osm_roads   | table  | user    |
| public   | planet_osm_ways    | table  | user    |
| public   | spatial_ref_sys    | table  | user    |
+----------+--------------------+--------+---------+

I do not like that the table names start with planet. Is it possible to define how the tables should be named?


texts: Running ftools multiparts to singleparts from PyQGIS

texts:
Is it possible to run the multiparts to singleparts function from the console with the processing_runalg function?
Or is there another way?


texts: How can I grey all outside a country border in Mapbox GL JS?

texts: I have klokantech docker container for openmaptiles. Now I want take the KML of the country border and grey all outside the country.


texts: Styling the map tips pop up box with css - background-color not working

texts:
I'm trying to find some extensive documentation on styling the map tips box when hovering over features. I'm using CSS school to do basic styling but some (like background) don't seem to do anything for me.

I'm using the following:

<style>
    h1 {color:red; font-size: 20px; text-align: center; width: 300px;}
    h3 {color:black;font-size: 12px; font-style: bold; font-family: verdana; border-width: 2px; border-left-style: solid;}
    p {color:black; font-size: 12px; font-family: verdana;}
    body {background-color: grey;}
</style>

So all of that works, except for the background. Among correcting my mistake, is there any QGIS specific css documentation other than https://nathanw.net/2012/08/05/html-map-tips-in-qgis/ which is not at lot to go off.


texts: Export GeoJSON to MSSQL with UTF-8 encoding

texts:
I'm testing this command to migrate GeoJSON data to MSSQL table:

ogr2ogr -f MSSQLSpatial "MSSQL:server=.\SQLEXPRESS;database=geodb;trusted_connection=yes" test.json -nln test

It works beautifully, but if json data contains other characters then ASCII I get garbage in MSSQL table for that field:

Json file is UTF-8 encoded without BOM header.

I tried appending -lco ENCODING=UTF-8 parameter to ogr2ogr command, but that didn't help.

So I wonder if anyone knows how to set sqlserver and ogr2ogr to communicate Unicode?


texts: Converting .mxd files to be compatible with earlier version?

texts:
I am running ArcMap 10.1. I need to open a .mxd that saved using 10.3, however, I do not have access to 10.3. I tried to execute a python tool that was created using Esri's How To: Save map documents in ArcGIS 10.x to previous version of ArcGIS in batch.

The error that I got was:

  Traceback (most recent call last): File
  "D:\GIS_data\scripts\mxd_rollback.py", line 14, in
  mxd.saveACopy(outloc + os.sep + base, version) File "c:\program files
  (x86)\arcgis\desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
  return fn(*args, **kw) File "c:\program files
  (x86)\arcgis\desktop10.1\arcpy\arcpy_mapping.py", line 826, in
  saveACopy self._arc_object.saveACopy(file_name, version)
  AttributeError: MapDocObject: Unable to save. Check to make sure you
  have write access to the specified file and that there is enough space
  on the storage device to hold your document.

  Failed to execute (ConvertMXDs).

Following is a copy of the script that I executed:

import arcpy, sys, os, string

mxdList = string.split(arcpy.GetParameterAsText(0), ";")
outloc = arcpy.GetParameterAsText(1)
version = arcpy.GetParameterAsText(2)

suffix = "_"+ version.replace(".", "")

for item in mxdList:
    item = item.strip('\'')
    mxd = arcpy.mapping.MapDocument(item)
    base = os.path.basename(item)
    base = os.path.splitext(base)[0] + suffix + os.path.splitext(base)[1]
    mxd.saveACopy(outloc + os.sep + base, version)
    arcpy.AddMessage(os.path.basename(item) + " has been converted")


texts: How to programmatically fill a datetime field with datetime.datetime.now() function

texts:
I have been trying to programmatically fill a datetime field using values coming from datetime.datetime.now() function.

#I create the temp shp using a URI
URI = 'MultiPolygon?crs=epsg:3857&field=id:integer&field=x:integer&field=y:integer&field=z:integer&field=bx:integer&field=by:integer&field=bz:integer&field=rating:integer&field=timestamp:datetime&index=yes'
# temp shp
mem_layer = QgsVectorLayer(URI, 'temp', 'memory')

# when I try to fill out the timestamp field with a variable now
now = datetime.datetime.now()
# I get an empty or zero value field

It looks like I have to cast the variable first to a string to make it work but in that case I lose the capability to make calculation over the field (I do need to keep track of time to evaluate the "duration" of a specific action triggered by the user)

I tried with the following datatypes with no luck:

date

datetime

real

double

Which datatype should I use when creating the shp in order for it to be filled with that specific value?


texts: Batch processing Surface volume calculation in ArcMap?

texts:
I'm trying to collect the 2D area, 3D area and volume of a given basin at 0.5 meter increment in plane height, but since I want to graph the relationship between plane height and the areas and volume, i want to collect many many points.

Doing this one by one is tedious. I'm wondering whether there is a python code that allows me to process a whole bunch of plane heights at once and congregate results in one .txt file or excel file.

currently, the code I'm using is

import arcpy
text="c:\\document\riverside\text.txt"
arcpy.SurfaceVolume_3d("watersheddem",text,"BELOW",250,1,1)

but this code only allows me to read result on the side bar, when i actually check the text.txt file it is empty even when my arcpy.surfacevolume_3d was ran successfully.


texts: What does this error means less than 4 coordinates in ## polygon

texts:
I'm using the bkde2D {KernSmooth} function on a dataset that has no repetitions of coordinates. while trying to present the data on a map I get this error "less than 4 coordinates in ## polygon".
.In a different dataset, there is no problem so I guess it's related to the current data. How can I overcome this problem?


texts: Setting spatial reference for MapInfo tab-file

texts: How can I change spatial reference for tab-file in MapInfo Professional?


texts: Pyqgis 2.2 : How to allow editing on layer and features

texts:
I'm writing  a pyqgis plugin with qgis 2.2
I succeeded in accessing and displaying vector data (geometries & attributes).

Now I would like to allow editing for all my layers and all features with layer’s in-memory editing buffer in order to be able in a 2nd time to commit/rollback the editions :

-> modifying features (geom & attributes)
-> adding features (geom & attributes)
-> deleting features (geom & attributes)

I read the qgis cookbook part called "Modifying Vector Layers with an Editing Buffer". But there is not examples showing the structure of the code. I would like to get, if possible the method I should use? Is it possible to get the structure of the code please? Which classes and which functions should i use ?

Thanks.

I'm writing a pyqgis 2.2 plugin with qgis 2.2 I succeeded in accessing and displaying vector data (geometries & attributes). I succeed also to create a 'beforecomitechanges' signal ( emitted, before changes are commited to the data provider). Thus a message has appeared into qgis logwindow ("comit change" label) after editing operations and before commiting toward the dp). Please read following piece of code allowing to do that:

def beforeLayerCommitChange(self):

QgsMessageLog.logMessage( "commit change" )

QApplication.beep()

layer.beforeCommitChanges.connect(self.beforeLayerCommitChange)

Now I would like to generate others signals. i would like generate the displaying of information messages (number of features that have been added,modifyied,deleted, the id of the features which have been modified, the name of the layer..).

COuld you throw light for me ? THanks.


texts:
I've got an XY dataset that I'd like to convert to a raster, based on the ETrel_whol field. The ETrel_whole field includes NaNs and ranges from 0-1. I would like to convert this to a raster, with a color stretched from 0-1, and NaN mapped as no color. The columns of the input data are shown here:

I convert to a raster using the Feature to Raster tool in the Conversion toolbox:

Note that I've also tried the Point to Raster, with the same outcome. It looks like ArcMap 10 assigns 102 "unique values" to the ETrel_whol field. 1 for NaN, and then one for each value from 0-1 at 0.01 intervals, based on my input data. The output attribute table is shown here:

The issue is that the Value field appears to be assigned randomly to my ETrel_whol field, rather than sequentially, and the output raster has the colormap based on the Value field, rather than the ETrel_whol field. If I switch my Symbology to Unique Values, I can define that the value field should be ETrel_whol, and colors are assigned to each value from 0-1, as shown here:

However, if I want to set the symbology to Stretched, it does not give me a choice of the value field, and simply stretches the Values from 1-102, rather than ETrel_whol from 0-1, as shown here:

How can I map the colors of this raster to the ETrel_whol field, rather than the Value field?

Thanks for your help!

-sam


texts: how do I join spatially polygons to polygons, creating new rows as required?

texts:
I have two polygons layers - counties and urban-mega-messes (UMM). I have a shapefile for each. I would like to spatially join the county attribute information to the UMM layer and create a new joined file, answering: what counties are the UMM's located in?

In many cases, the UMM's are in multiple counties, in which case I would like additional rows to be added to the joined shapefile for each county in which they are contained.

UMM 1     COUNTY 1
UMM 1     COUNTY 2
etc

For example, the UMM "New York City" is in multiple counties.

I have been trying Data Management Tools -> Join attributes by location, but only 1 county record is appended to each UMM in the joined table - no additional rows are added.

I've searched online but found no answer. Is there a plugin that would do the trick? Or another tool in QGIS?


texts: Changing text size or line height of QGIS attribute table?

texts: How do I change the text size in the attribute table or the height of the rows in QGIS?


texts: Combined area of all features within features of another layer

texts:
I have a grid shapefile (basically, a half-degree grid of the world, where each cell has a unique id number), and a set of files with features. I want to record for each grid cell, all the areas and/or values of one field of the second set of files that fall within the gridcell. The "fall within" is whether the centroid is within the cell. My code looks like this (Python+OGR, I can't get shapely installed on the work system)

    # grid_layer is my grid layer, it has been opened elsewhere
    # outut_result is a dict indexed by grid ID#
    feat_shape = ogr.Open ( shapefile )
    feat_l = feat_shape.GetLayer ( 0 )
    extent = feat_l.GetExtent()
    # Apply the extent filter to the grid layer
    grid_layer.SetSpatialFilterRect ( *extent )
    # Now loop over each grid cell within the feat region
    for grid in grid_layer:
        grid_geom = grid.GetGeometryRef()
        grid_no = grid.GetFieldAsInteger ('ID')
        # Loop over features within this grid id
        # Start by applying a spatial filter
        feat_l.SetSpatialFilterRect ( *grid_geom.GetEnvelope())
        for feat in feat_l:
            # Get the feat's geometry and calculate its centroid
            feat_geom = feat.GetGeometryRef()
            centroid = feat_geom.Centroid()
                # Is the centroid within the grid_geom?
                if centroid.Within ( grid_geom ):
                    value = feat.GetFieldAsInteger ( 'field_of_interest' )
                    # Or feature's area, or whatever
                    # Store it in our output dict
                    output_result.setdefault(grid_no, []).append( value )

            del centroid
            del feat_geom

        del grid_geom

At the end of this process, I am left with output_result and my information: output_result[29345] results in a list with all the values of field_of_interest, which I then can process further with e.g. numpy or scipy.

My question is that this processing takes forever, as there are a lot of shapefiles. Then writing it out is rather messy (doing it with a pickle at the moment, but that needs even more code afterwards). Am I missing some really obvious "point-and-clicky" way of solving this problem?


texts: How to open G.L.O.B.E. DEM files from NOAA

texts: I need to open a free GLOBE DEM. I have downloaded some tiles from: https://www.ngdc.noaa.gov/mgg/topo/gltiles.html, but I don't know which software should I use in order to open those files. Is ArcMap able to view those files ?


texts: How to update OSM planet diffs using flat nodes file?

texts:
I'm importing the OSM planet to my PostgreSQL using the flat file option.

Now I have my database and my flat nodes file, I want to apply some diffs files.

After reading this :

  Lastly, since you have to reimport for --flat-nodes, you should make sure to download a new Europe extract, and make sure to download PBF, not bzipped XML.

I understood I need to download the diff file and ALL planet file again. Is it correct? By importing without the flat file option and keeping all temp tables in database, I just need the diff file to do the update.

Please explain the correct process to apply diff updates using the flat nodes file.


texts: Correct SQL expression for SelectbyAttribute Tool

texts:
I'm trying to apply the sql expression described here https://geonet.esri.com/thread/76028, that involves querying out odd values in a field. The specific sql expression is Round(NumberField / 2, 0) <> NumberField / 2.

I've written this in python:

import arcpy

plantFile = r"U:\Users\K\Plants.shp"
field = "UID"

expression = Round("UID" / 2, 0) <> "UID" / 2
arcpy.SelectLayerByAttribute_management(plantFile, "NEW_SELECTION", expression)

and it returns this error: NameError: name 'Round' is not defined

The expression works when applied in ArcGIS, but not in arcpy. Any hints?


texts: ESRI Javascript v API 3.3: Format Info Window Content

texts:
I'm attempting to query a parcel layer and display the information in a popup infoWindow when I click on a parcel, and I'm having some issues dealing with different datatypes.

var content = "Address: ${ADDR}" + "
Owner Name: ${OWNER}" + "
Parcel ID: ${APN}" + "
City: ${CITY}" + "
Acres: ${TOTAL_ACRES:NumberFormat}" + "
County Assessor Site"
var popUpTemplate = new esri.InfoTemplate("Parcel", content);

The ADDR, OWNER, APN and CITY datatypes are all string so they appear in the infoWindow without issue, but I'm having trouble getting TOTAL_ACRES and COUNTY_LIN to appear.

TOTAL_ACRES is a double, and i thought that using the dojo.number.format method as outlined on the esri website would help, but it doesn't seem to do anything.

COUNTY_LIN is a url, and has a data type of text, but it isn't appearing and I'm not exactly sure why.
I've included screenschot of the relevant data info from inside ArcGIS, hopefully they are helpful. Any help is appreciated! I'm new at this and very stuck as to what I'm doing wrong.


texts: I wrote a simple script in VB. But it failed to run in field calculator. I can't figure out what is wrong with it. Below is a screenshot.


texts: is there an API-level error handler in the ArcGIS Server JS API?

texts:
I'm seeing an error message in FireBug:

Error: Internal kml parser error.http://serverapi.arcgisonline.com/js.../?v=3.2compact
Line 34

I'm not worried about fixing this error for the moment - I'd like to know how can I trap the error, so I can report it to the user elegantly? Eg, if they're not using Firebug I'd still like them to know that there's a problem.

The map has an onLayersAddResult event, which contains an Error handler. Is there a similar map-level error handling event which would trap the above KML parser error (and any other errors which might be generated by the API)?

It seems like the map object needs a generic onError() event, but I can't see anything suitable.

Thanks


texts:
Why do the QuickMapServices plugin generates distorted basemap?
I tried different sources but the result was odd every time.
I'm using QGIS 2.16.2 and QuickMapServices 0.17.8.


texts: Recorded track capacity of Garmin eTrex n0 devices

texts:
I am interested in purchasing a Garmin eTrex 10 or 20 primarily for geotagging the photos I take on holiday, although I suppose the other features might come in handy as well in the future. When using the device for geotagging I would (like to) record my movements during a day and save them to a .gpx file and copy the .gpx file to my laptop at the end of the day when I process my photos, of course I would clear the storage of the Garmin then too. I think a track point interval of 5 m or 5 s would suffice for my use. Unfortunately I am a bit confused about the capacity of the eTrex n0 devices. According to the Garmin website a/the track log can consist of 10 000 track points on both the eTrex 10 and 20(x). However it isn't clear to me what this means. Does this mean every track can consist of not more than 10 000 track points? Or is this a total?
Several people on several forums have claimed that this is the total amount of track points for all tracks the device(s) can store.

Edit:
According to Outdoor Gear Lab - The Best Handheld GPS Review the 10 000 is in fact the amount of points per track.


texts: Populating ESRI GeoPortal programmatically?

texts:
I'm trying to find either some documentation or an example of how to manually populate the GeoPortal metadata datastore.

I'm using a MySQL backend.  Given a ISO XML file I can manually harvest it but this is tedious and doesn't fit with our model. I can easily write a script to populate entries in the SQL database but it appears that there is more to it than just adding entries to SQL, there is a part that goes into the Lucene engine as well.

Does anyone have any experience or can provide a simple example, how-to or preferably some documentation on populating the datastore?

The goal is to have GeoPortal as an external catalog of locally registered datasets, the trouble is registering the datasets with the GeoPortal catalog in a preemptive manner rather than a reactive (having an administrator manually "harvest" or ingest resources).


texts: Create centroids from point layer based on unique field names

texts:
I have a point layer in QGIS, which I want to filter by attribute values of a specific attribute column. For each unique attribute value, a centroid should be calculated from the respective features. This results in one centroid for each set of features containing a unique attribute value.

Example:
In the attribute column 'names', I have 10 features named 'aaa', 20 features named 'bbb', and so on. Now I want to create a centroid of the respective 10 'aaa' features, another centroid of the respective 20 'bbb' features, etc.

Anyone has some hints on how to do this?


texts:
I have a list of points and a list of polygons both in a separate shape files. I am trying to write a ArcPy script to loop through the points, and then the polygons respectively and tell me when a point is within a polygon. I have tried both the contains and touches methods on polygons but haven't had any luck. Any help would be appreciated. So far I have the following:

    pointFeatureClass = "C:\\Data\\POINTS2.dbf"
    pointFieldList = arcpy.ListFields(pointFeatureClass)
    for field in pointFieldList:
        curField = field.name
        arcpy.AddMessage(curField)

    ptRows = arcpy.SearchCursor("C:\\Data\\POINTS2.dbf", "", "", "CID;Shape","CID;Shape")
    for ptRow in ptRows:
        #Now loop through poly's testing for intersect
        shpRows = arcpy.SearchCursor("C:\\Data\\counties.shp", "", "", "NAME;Shape","NAME;Shape")
        for shpRow in shpRows:
            #arcpy.AddMessage("Checking for point "+ str(ptRow.CID) +" in: " + shpRow.NAME)
            if shpRow.Shape.touches(ptRow.Shape):
                arcpy.AddMessage("Got it: " + shpRow.NAME)


texts: Carto Builder, show attribute readout synced to time data in animated map?

texts:
I'm working on a donated webmap animating time data, the progress of small electric boat across the ocean. The time field for animation is a time in seconds counter from launch attribute. The data table also has fields for proper date (text as mo,day,year) and boat telemetry for solar voltage, pitch, and roll.

Is it possible to have readouts in the map for these other attributes, updated by the specific time point in the animated data?


texts: Contents of Downloadable Data from ESRI website?

texts: I'm looking for the simple grey basemap to use in a map I'm working on, but want to download it.  It's one of the options with their dynamic data.  ESRI has no real order on their download page.  Anyone with a list of contents inside these data downloads? I figure the newer ones just have more updated data.


texts: Is there any way to map match static points to roads using PostGIS

texts:
I'm using PostGIS and Leaflet for positioning elements on map. Because of GNSS (GPS) error, points are not always on road. Is there any way to eliminate that error and position element to road?
I don't have roads in database (for now).

Elements are points with gps coordinates and some attributes. Is there any way to map match elements?


texts:
I am trying to define a custom CRS based on Equal Earth Projection (As described here) in QGIS (Settings >> Custom Projections >>  + [Add new CRS] >> Entering Proj String [+proj=eqearth +datum=WGS84 +wktext] >> OK).
The validation of the Proj String goes fine; however, after pressing OK it gives the following error message.

Cannot save' * Generated CRS (+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 + datum=WGS84 + units= m + no_defs)'
— this Proj string definition is identical to ESRI:54032.
Try changing the CRS definition to a WKTformat instead.

I thought that probably ESRI:54032 and Equal Earth maybe the same thing. But when I change my map projection to ESRI:5432, it does not looks like Equal Earth Projections.

Are the Equal Earth and ESRI:54032 the same?
If not, how can I define the Equal Earth projection in QGIS?

I am using QGIS 3.12 (64 bit) on Windows 10. Detailed specification about QGIS and related libraries is here.

A similar question already exists but neither the question nor the answers address the problem I have described.


texts: Adding new reclassified column to attribute table in QGIS using "CASE", "WHEN" "THEN"

texts:
I am stuck by adding new column (kg/ha) based on "Value" table with using following. This is an NDVI layer which I want to convert to prescription map for variable rate fertilization. I tried following:
case
when"VALUE" <= 0 then 20
when"VALUE"= 0-0.05 then 30
when"VALUE"= 0.05-0.10 then 45
when"VALUE"= 0.10-0.15 then 130
when"VALUE"= 0.15-0.20 then 125
when"VALUE"= 0.20-0.25 then 120
when"VALUE"= 0.25-0.30 then 110
when"VALUE"= 0.30-0.35 then 95
when"VALUE"= 0.35-0.40 then 85
when"VALUE"= 0.40-0.45 then 80
when"VALUE"= 0.45-0.50 then 75
when"VALUE"= 0.50-0.55 then 70
when"VALUE"= 0.55-0.60 then 65
when"VALUE"= 0.60-0.65 then 60
when"VALUE"= 0.65-0.70 then 55
when"VALUE"= 0.70-0.75 then 50
when"VALUE"= 0.75-0.80 then 45
when"VALUE"= 0.80-0.85 then 40
when"VALUE"= 0.85-0.90 then 35
when"VALUE"= 0.95-0.95 then 30
when"VALUE">=0.95 then 25
end

but the result is only "NULL". Can someone tell me what I doing wrong?


texts: Merging Fields into new field using python and eliminating space in before and after the field using field calculator

texts:
Am interested in combining ClASS1 and SUBCLASS1 to form the one combined class.  I wrote the expression as follows "() ()" format( !CLASS1!replace(" "","" ),( !SUBCLASS1!replace(" "","" ))as shown in the Image[![field to be combined and expression] [![enter image description here]

but after executing I get the following error message


texts: ogr2ogr, force DATE cast on a SQL request

texts:
I am trying to complete column null values on a .shp file with one default value. My column has a DATE type.

Files schema : "COLUMN_NAME TYPE(SIZE)"

Input (with null values) :

CLICHE VARCHAR(80) | DATE DATE(10) | HEURE_TU VARCHAR(80)

Output (with defautl value) :

CLICHE VARCHAR(80) | DATE_PVA DATE(10) | HEURE_TU VARCHAR(80)

I tried:

ogr2ogr -dialect sqlite -sql "SELECT geometry, CLICHE, CASE WHEN DATE IS NULL THEN '2012-01-01' ELSE DATE END as DATE_PVA, HEURE_TU FROM 'my_table'" outputfile.shp inputfile.shp

But I get a file with a VARCHAR(80) typed DATE_PVA column.

So I tried:

ogr2ogr -dialect sqlite -sql "SELECT geometry, CLICHE, CAST((CASE WHEN DATE IS NULL THEN '2012-01-01' ELSE DATE END) as DATE) as DATE_PVA, HEURE_TU FROM 'my_table'" outputfile.shp inputfile.shp

But I get a output file with a INT(9) type DATE_PVA column.

I am pretty sure that the 'aaaa-mm-dd' format is the good one, however I have also tried to modify date format '01-01-2012' and '01/01/2012'. No better results...

So I do not know what I can do. I need a command line solution to be able to run it in a script. However, I can try an other solution than ogr2ogr if needed.


texts:
I am writing a script to create thousands of maps and I am trying to create a table that updates for every map. Basically an updating table of the information specific to each map. I have been trying to link a an excel table to the layout view but arc does not seem to have that functionality. So what I am doing instead is creating a manual table in the layout view. I have 27 values for each map that I have been storing in a dictionary (using search cursors to extract the info and store into dictionary).
I cannot copy and paste my table into here but it essentially looks this:
Receptors               AEGbuffer    CONbuffer    ERGbuffer
buffer                   1            10           19
schools                  2            11           20
childcares               3            12           21
hospitals                4            13           22
nursinghomes             5            14           23
critical                 6            15           24
sw                       7            16           25
streams                  8            17           26
respop                   9            18           27

for k,v in sorted(table.items()):
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if elm.text == "1":
            if k == 'AEGbuffer': elm.text= v
        if elm.text == "2":
            if k == 'AEGschools': elm.text = v
        if elm.text == "3":
            if k == 'AEGchildcares': elm.text = v
        if elm.text == "4":
            if k == 'AEGhospitals': elm.text = v
        if elm.text == "5":
            if k == 'AEGnursinghomes': elm.text = v
        if elm.text == "6":
            if k == 'AEGcritical': elm.text = v
        if elm.text == "7":
            if k == 'AEGsw': elm.text = v
        if elm.text == "8":
            if k == 'AEGstreams': elm.text = v
        if elm.text == "9":
            if k == 'AEGrespop': elm.text = v

this is just for the first row, the rest of the code is the same process for the pasting of the dictionary values.
where I am stuck: I do not know how to properly reset the newly changed values to their original 1-27 values so the next map can be easily produced. is there a tool or trick to this?
what I have tried: after I export the map I reopen the text elements and reverse the process like..
for k,v in sorted(table.items()):
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
         if k == 'AEGbuffer':
             if elm.text == v: elm.text = "1"
             and so on......

but this approach hasn't worked great. a few of the values are wrong when there are the same values a few times on the table.


texts: Zoom to current features in SQL layer in OpenLayers 4.6?

texts:
I've created a SQL layer, I'm passing the values for that query from client, I want my map to be able to zoom to whatever the layer that comes after query. I tried GetCapabilities, but it zooms to the extent of whole layer, which is not quite helpful.

I've many layers on GeoServer which, when I click on checkbox, displayed on the map.

This is what I'm using now in my program

var url = 'http://localhost:8080/geoserver/wms?;
request=GetCapabilities&service=WMS&version=1.1.0';
var parser = new ol.format.WMSCapabilities();

$.ajax(url).then(function(response)
{
    var result = parser.read(response);

    var Layers = result.Capability.Layer.Layer;
    for(var i=0; i<Layers.length;i++){
     var layerobj = Layers[i];

      if (layerobj.Name == 'SPA:'+value.id) {

              var ex = layerobj.BoundingBox[0].extent;
              //ext = map.getView().calculateExtent(map.getSize());

              var X = ex  [0] + (ex  [2]-ex  [0])/2;
              var Y = ex  [1] + (ex  [3]-ex  [1])/2;
      }
    }
    map.getView().setCenter([X,Y])
    map.getView().setZoom(19);

Using this code, I'm able to set zoom level and center of map according to complete layer, but that is not what I want, I want to set zoom level and center according to the features that I get after query.


texts: Importing/downloading symbology from ArcGIS Online using ArcPy?

texts:
I am currently attempting to create a script to download shapefiles from ArcGIS Online to my local drive before adding it to an aprx project and thereafter uploading it to a ArcGIS portal/enterprise.

All the above I have managed by downloading them as shapefiles -> convert .shp to feature class (in gdb) and then uploading it to the aprx.

However, when I run the script the symbology do not tag along. I have tried finding methods of downloading lyrx files from AGOL without any luck.

Do any of you know of a way to download the symbology from AGOL and merge them with the shapefile or feature class?

Here's the most relevant functions in the code I've used:

def agol_login(AGOL_url,AGOL_user,AGOL_pw):
    # gis_q = AGOL_url + ', ' + AGOL_user + ', ' + AGOL_pw
    print('AGOL login')
    gis = GIS(AGOL_url,AGOL_user,AGOL_pw)
    return gis

def createGDB():
    if folder_exist(path to gdb):
        gdb = 'testpath'
    else:
        gdb = arcpy.CreateFileGDB_management(path to folder, gdb)

def folder_exist(path):
    if arcpy.Exists(path):
        print(path + " Exist")
        return True
    else:
        print(path + " does not exist")
        return False

def downloadUserItem2(item, export_location, downloadFormat):
    # type: (object, object) -> object
    try:
        #Download
        name = item.title.replace(" ", "_")
        result = item.export(name.format(item.type), downloadFormat)
        result.download(export_location)
        result.delete()

        #Unzip
        zip_path = export_location + "\\" + name + ".zip"
        zip_file = ZipFile(zip_path)
        zip_file.extractall(export_location)

    except Exception as e:
        print(e)

def SHPtoGDB():
    # Set the workspace for ListFeatureClasses
    arcpy.env.workspace = "path to folder"
    # Use the ListFeatureClasses function to return a list of shapefiles.
    featureclasses = arcpy.ListFeatureClasses()
    # Copy shapefiles to a file geodatabase
    for fc in featureclasses:
        arcpy.CopyFeatures_management(
            fc, os.path.join("path to gdb",
                             os.path.splitext(fc)[0]))

def AddFCtoAPRX():
    arcpy.env.workspace = "path to gdb"
    datasets = arcpy.ListDatasets(feature_type='feature')
    datasets = [''] + datasets if datasets is not None else []
    rootPath_aprx = 'path to aprx folder'
    aprx = arcpy.mp.ArcGISProject(rootPath_aprx + 'TestWMS.aprx')
    aprxMap = aprx.listMaps("Map")[0]
    for ds in datasets:
        for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
            path = os.path.join(arcpy.env.workspace, ds, fc)
            print ("List of feature classes:")
            print(path)
            try:
                test = aprxMap.addDataFromPath(path)
                aprxMap.addLayer(test)
                print("Feature Classes added to Arcgis Pro..")
            except:
                print("gah.")
    aprx.save()
    print("APRX path is:")
    print(aprx.filePath)

def UploadtoPortal():
    rootPath_aprx = 'C:\\Users\\ESSE\\Documents\\ArcGIS\\Projects\\TestWMS\\'
    aprx = arcpy.mp.ArcGISProject(rootPath_aprx + 'TestWMS.aprx')
    m = aprx.listMaps('Map')[0]
    arcpy.mp.CreateWebLayerSDDraft(m, rootPath_aprx + 'TestWMS.sddraft', 'Map', 'MY_HOSTED_SERVICES', 'FEATURE_ACCESS')
    arcpy.StageService_server(rootPath_aprx + 'TestWMS.sddraft', rootPath_aprx + 'TestWMS.sd')
    print("Uploading Service Definition to Portal...")
    arcpy.UploadServiceDefinition_server(rootPath_aprx + 'TestWMS.sd', 'My Hosted Services')


texts:
I want to mosaic 2 rasters "rasa" and "rasb" using gdal_merge.py . While running the code, it produces an error " No input files selected". Can somebody help me to find the error in the code?

    import subprocess
    filea = 'C:/Users/claudio/workspace/test/test1/rasa.tif'
    fileb = 'C:/Users/claudio/workspace/test/test1/rasb.tif'
    output = 'C:/Users/claudio/workspace/test/test1/output.tif'
    subprocess.call(['gdal_merge','-o',output,filea,fileb],shell=True)

when i run the above code, the output i am getting is:


texts: Problem with adding placemark from kml file with specific id

texts:
I have a problem with adding a vector layer from kml file. I am trying to make XMLHttpRequest and add only a placemark with specific id. Nothing is displaying on the map. Request is successful and variable data contains proper snippet of kml file. I'm pretty sure I'm making some noob mistake. Here is my code:

map = new OpenLayers.Map("mapdiv");
map.addLayer(new OpenLayers.Layer.OSM());

var idOfPlacemark = "test";
var req = new XMLHttpRequest();
req.open('GET', 'boundaries.kml', false);
req.send(null);
if(req.status == 200) dump(req.responseText);

var data = req.responseXML.querySelector("[id=" + idOfPlacemark + "]");
var format = new OpenLayers.Format.KML;
test_layer = new OpenLayers.Layer.Vector("test");
test_placemark = format.read(data.outerHTML);

test_layer.addFeatures(test_placemark);
map.addLayer(test_layer);

var layer_switcher= new OpenLayers.Control.LayerSwitcher({});
map.addControl(layer_switcher);

var lonLat = new OpenLayers.LonLat(19.274483, 50.203275)
    .transform(
        new OpenLayers.Projection("EPSG:4326"), // transform from WGS 1984
        map.getProjectionObject() // to Spherical Mercator Projection
    );
var zoom=13;
map.setCenter (lonLat, zoom);

And here is my kml file:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.0">
<Document>
<Placemark id="test">
    <name>Name of placemark</name>
    <description>This is descr</description>
    <Polygon>
    <extrude>1</extrude>
    <altitudeMode>relativeToGround</altitudeMode>
    <outerBoundaryIs>
        <LinearRing>
            <coordinates> 19.212462, 50.200544
                        19.21287, 50.200399
                        19.21346, 50.19998
                        19.213857, 50.199603
                        19.214146, 50.199548
                        19.214715, 50.199541
                        19.215766, 50.199761
                        19.216743, 50.199987
                        19.217311, 50.200083
                        19.218041, 50.200028
                        19.21876, 50.199871
                        19.219178, 50.19961
                        19.219446, 50.199376
                        19.219253, 50.199046
                        19.218331, 50.198318
                        19.217741, 50.197845
                        19.217333, 50.197508
                        19.217086, 50.196986
                        19.216421, 50.196279
                        19.215606, 50.196004
                        19.214737, 50.196024
                        19.213996, 50.196182
                        19.213331, 50.196423
                        19.212065, 50.197055
                        19.211089, 50.1977
                        19.210874, 50.198298
                        19.210896, 50.199005
                        19.211217, 50.199747
                        19.211647, 50.200228
                        19.212054, 50.200592
                        19.212301, 50.200598 </coordinates>
            </LinearRing>
    </outerBoundaryIs>
</Polygon>


texts: Removing pixel values below 0 in Digital Elevation Model in ArcGIS for Desktop?

texts:
I made a DEM layer from a .las file (LiDAR).

How can I remove the black area which contains negative values? Is it  possible to set these values equal zero or null?


texts: Alternative to Terrasync for Trimble Windows Mobile unit

texts:
I have a number of Trimble Nomad 900 series devices and am currently using TerraSync.

TerraSync does the job, but are there any other alternatives out there that would be even simpler to use?

Basically I'm looking for something that's quick and easy to record features and save to a non-post processed file that can be named.  Preferably just a text file with the "average" or best coordinate reading.  The TerraSync workflow (i.e. having to use GPS Pathfinder Office) is overkill for our uses and takes too much time.


texts: Dynamically change omnivore file with global variable?

texts:
I'm trying to create a number of buttons that will switch out the CSV files I'm filtering rather than having to create separate layers for each file, as they all have the same structure. This seems plausible as I can define my CSV file as an initial variable and successfully call it with `window.dataMV = omnivore.csv(myvariable)", but whenever I toggle the button to change the value to File_B, this doesn't seem to affect the map.

Even refreshing the layer by removing window.dataMV and re-adding it with another button, still seems to show File_A, the initial value the global variable was given. Is there an easy workaround that allows Omnivore files to be changed dynamically with this basic Javascript click function?

   myvariable = 'File_A.csv';

$('#button').click(function() {

if(!toggle) {
    myvariable = 'File_B.csv';
    console.log(myvariable);

  } else {

    myvariable = 'File_A.csv';

  }
  toggle = !toggle;
});

console.log(myvariable);

window.dataMV = omnivore.csv(myvariable)
    .on('ready', function(layer) {
    markerArray = [];
    this.eachLayer(function(marker) {
        if (marker.toGeoJSON().properties.Category == "SAMPLE 1") {
            marker.setIcon(L.icon({
                iconSize: [18, 18],
                iconAnchor: [9, 27],
                popupAnchor: [1, -24],
                iconUrl: 'c1.png'
            }));

        else {
            window.dataMV.removeLayer(marker)}

      return marker.bindPopup(marker.toGeoJSON().properties.Category + ', ' +
                marker.toGeoJSON().properties.Location);
    });
//.addTo(map);
});


texts:
It has been a while since I have done anything but basic work with GIS.

I have a map with a utility system, comprised of points and lines. I wish to create a map with only a select few of these lines.

Is there any easy way to copy the vectors and nodes I want, and copy them to a new layer or shapefile, or to mass delete the one I do not want?


texts: Batch add multiple columns to a shapefile

texts:
I am working on a process and during that process for each project I need to add 16 new fields as text strings to my shapefile. T_1, R_1, S_1, QS_1,... these 4 fields out to "_4" .

My scripting knowledge is not quite up to par but learning.


texts: How a draw a circle (on the surface of earth) on a mercator projection?

texts:
I would like to use R and leaflet to draw a circle (a real circle on the surface of the earth) on a mercator projection.

In the following code, I have the center of the circle (with lng=0 and lat=60), and I would like to have a parametric equation to draw the circle. I started with the equation below, and I can draw it on a leaflet map:

centre=c(0,60)
t=seq(0,2*pi,0.005)
r=9
xx=cos(t)*r+centre[1]
yy=sin(t)*r+centre[2]

Then I can draw the form on a leaflet map:

leaflet() %>%
  addTiles() %>%
  addPolygons(
    lng=xx, lat=yy,fillColor = "transparent"
  )%>%
  addCircleMarkers(lng=centre[1], lat=centre[2],radius=1,color="red")

The code above gives this kind of formes on the earth:
http://www.mylovedone.com/image/solstice/sum12/tissotFig3b.jpg

I did't make the image, but I think they are the same. And what I want to draw is this: http://imgur.com/wLptl

How should I modify the equation of the circle so that I can obtain a real circle.

Is this a Tissot's indicatrix? It seems that a Tissot's indicatrix is to render a circle on the sphere on a map. However, the Tissot's indicatrix on a mercator projection is a circle! I can't understand how it is possible. Since if you draw a circle on the north pole, there is no way it is a circle on mercator projection.

EDIT: I found this web site http://www.movable-type.co.uk/scripts/latlong.html it shows the equation of destinations, givens one point, the direction and the distance. So it is then possible to draw the cercle!


texts:
I have a raster band saved as a .jp2 which I downloaded from Sentinel2-L2A. I can load it using rasterIO:

with rasterio.open(path_to_jp2) as src:
    img = src.read(1)

This works fine, however, it takes a long time to load since the tile is big. I'm wondering if there's any way to load only a specific section of the image, let's say a small field stored as a geopandas dataframe which I extracted from a .shp file?

field = geopandas.read_file(path_to_shp)

I was thinking about converting a .jp2 into a cloud optimised geoTiff but I don't know if that would be the most efficient solution. Any suggestions?


texts:
I installed Apache2 and GeoServer, then set up the proxy.cgi. I imported shapefile based multipolygon file into GeoServer. I can access this via the GeoServer OpenLayers viewer but when I want to do the same in an HTML file I can't do it, because the GetFeatureInfo doesn't want to work properly.

I added to the file:

OpenLayers.ProxyHost = "proxy.cgi?url=";
map.events.register('click', map, function (e) {
                    document.getElementById('nodelist').innerHTML = "Loading... please wait...";
                    var params = {
                        REQUEST: "GetFeatureInfo",
                        EXCEPTIONS: "application/vnd.ogc.se_xml",
                        BBOX: map.getExtent().toBBOX(),
                        SERVICE: "WMS",
                        INFO_FORMAT: 'text/html',
                        QUERY_LAYERS: map.layers[0].params.LAYERS,
                        FEATURE_COUNT: 50,
                        Layers: 'csongrad_demo:csongmat',
                        WIDTH: map.size.w,
                        HEIGHT: map.size.h,
                        format: format,
                        styles: map.layers[0].params.STYLES,
                        srs: map.layers[0].params.SRS};

                    OpenLayers.loadURL("http://xxx.xxx.xxx.xxx:80/geoserver/csongrad_demo/wms", params, this, setHTML, setHTML);
                    OpenLayers.Event.stop(e);
                });

Then I tried this:

info = new OpenLayers.Control.WMSGetFeatureInfo({
                    url: 'http://xx.xx.xx.xx:80/geoserver/wms',
                    layerUrls: 'http://xx.xx.xx.xx:80/geoserver/csongrad_demo/wms',
                    title: 'Identify features by clicking',
                    queryVisible: true,
                    eventListeners: {
                        getfeatureinfo: function(event) {
                            map.addPopup(new OpenLayers.Popup.FramedCloud(
                                "chicken",
                                map.getLonLatFromPixel(event.xy),
                                null,
                                event.text,
                                null,
                                true
                            ));
                        }
                    }
                });
                info.activate();
                map.addControl(info);

None of them wanted to work, I always get this error "Loading... please wait..."


texts: Linearized map of rivers - what is this called?

texts: I am looking for a "linearized overview" map (as in the figure below) of the Brahmaputra river. I don't know the exact name of this type of map. Do you have any suggestion on what I should look for, and where?


texts: Combining multiple bands from different image collection into one collection

texts:
I did some calculations on variety of image collections and have generated products and saved them in variables, say a, b and c. Now I want to combine all these products (All one bands) into a multiband image collection.

What function should i use?
I tried with:

a.combine(b).combine(c)
a.merge(b).merge(c)
a.add(b).add(c)

But, these functions are not working.


texts: Global grid, which reference system to use?

texts: I have seen a couple global studies using WGS84 as geographic reference, with the problem that the area of raster cells is not constant. Which reference system would you use for a global analysis in order to have a constant cell size?


texts: Perform a spatial join in C# with ArcObjects

texts:
I feel like this has to be on the internet somewhere, but I haven't been able to find it.

I'm attempting to perform a spatial join on a polygon and point layer.  My goal is to push the polygon's ID to each point that it contains.  In PostGIS, it might look close to this:  SELECT a.*, b.* FROM a, b WHERE ST_CONTAINS(a, b).

I have:

A polygon IFeatureClass, from a geodatabase
A point IFeatureClass, from a geodatabase
The name of the polygon field

I would like to turn this into a point IFeatureClass with an extra field.

My question: How do I do this?  Bonus points for links to relevant documentation.


texts: looking to integrate R-trees into AlchemyDB

texts:
I am the author of AlchemyDB: http://code.google.com/p/alchemydatabase/

I would very much like to add geospatial indexing in the form of R-trees (most likely) and need some advice on whether or not there any C/C++ packages that have simple stripped down R-trees. I am a big fan of small codebases that are easy to read.

If anyone has any suggestions or knows of any experts that I could bug (I can pay for their time), I would be grateful


texts: More efficient processing of many XYZ tiles into a merged response

texts:
I've built a mobile mapping device for my Search & Rescue group using a Raspberry Pi (3b+) and its primary purpose is to serve a cache of XYZ tiles via its web server, over its WiFi network. The device is powered using a vehicle's 12V power outlet. Our SAR group uses a number of mapping applications that understand the XYZ tile format. We are regularly outside cell service and satellite internet is problematic, so the goal is to take data with us when we head out on a call.
As a bonus feature I've added the ability to export a georeferenced PDF of any area, with the PDFs content built from the relevant XYZ tiles using GDAL. This works OK but it's slow and I'm wondering if I can remove any of the more expensive steps from the transformation.
The current steps are:

identify all XYZ tiles at a given zoom level that intersect the caller's lat/long bounding box (this is fast)
call GDAL Transform for every XYZ tile, converting it to a .tif with its outputBounds set to the correct bounding box. This can mean several thousand calls to Transform and this part is very slow
call GDAL Warp to merge all Transformed .tif tiles into a single .tif (relying on the newly-added georeferencing data from step 2)
call GDAL Transform to convert the merged .tif to a PDF

I'm already using the Pi's "overlay file system", which I believe writes exclusively to memory / swap space - it's effectively a temporary file system where all changes are lost on reboot - and this made the export significantly faster. This also helps to mitigate the risk of storage corruption as the power supply is regularly interrupted.
It seems like there should be a way to remove step 2, where hundreds or thousands of tiles are individually georeferenced before being merged, but this would require some way of telling GDAL a) the overall bounds of the combined XYZ tiles and b) the order in which the XYZ tiles should be assembled, because without georeferencing we don't know where each tile belongs in the world.
Can anyone suggest a way to make this process more efficient?


texts: How to get the envelope of all the features in a feature class in ArcMap?

texts: I would like to get a single envelope of all the features in a feature class. Is there any way to do that in ArcMap?


texts: How to compute max and min coordinates for multiple polygons?

texts:
I would like to calculate the maximum, median and minimum latidudinal coordinates for a large set of polygons, using QGIS or R. So the highest and lowest points reached by each polygon across a latidudinal range. Central coords are easy to do, but converting the polygons to points, and using a distance matrix in QGIS crashes my computer. Is there a more effecient way to do this?

Its been asked in various guises here, but either 1. Unclear 2. outdated I think version wise or 3. crashes on my PC as I have about 2000 polygons.


texts: Spatial file attribute table shows null after joining CSV file

texts:
While I join a data table (saved as CSV) and added as both ways- 'delimited Text layer' and after that as 'new vector layer'.

After matching layer, the data shows null in the spatial attribute file. An earlier version used to join it without hassle.

After following some tips, CSVT seems to be a solution but generating a sidecar file


texts: Mapping Active Directory groups to GeoServer roles

texts:
I am trying to grant the GeoServer ADMIN role to members of an Active Directory group.  I have followed the documentation here and here.  My Active Directory users can successfully login, but the members of my Active Directory admin group do not receive ADMIN privileges in GeoServer.

Is there a problem with my configuration (below)?

I have tried many, many other filter values, without any success.  Even if you can't offer a solution, comments from anyone who has successfully used Active Directory groups in GeoServer, will let me know the functionality works as documented.

Here are my LDAP Authorization settings in GeoServer (values have been changed for public posting, but the general structure is accurate):

Use LDAP groups for authorization: checked
Bind user before searching for groups: checked
Group search base: OU=Groups,OU=Department,OU=Company Name,DC=domain,DC=com
Group search filter: (member={1}) where {1} is apparently the user's complete distinguished name
Group to use as ADMIN: GeoServer Admins
Group to use as GROUP_ADMIN: GeoServer Admins

When I use an LDAP utility (Softerra LDAP Administrator), the following settings will successfully return a user's groups:

Search DN: OU=Groups,OU=Department,OU=Company Name,DC=domain,DC=com
Filter: (member=CN=Firstname Lastname,OU=Users,OU=Department,OU=Company Name,DC=domain,DC=com)

(I have not set up an LDAP Role Service in GeoServer, but I don't think I need to if I simply want to map one AD group to the ADMIN role in GeoServer?)


texts:
I have a dataset of rasters (driver kea) whose has, let's say, 20-30 attributes in their RAT. I'm realizing how to have these attributes in shapefiles (without using proprietary software). I've read a lot of stuff in GDAL and Python, but I was not able to solve it. Since I can't convert my rasters to shapefiles via gdal_polygonize.py (It doesn't deal with RATs), I'm trying to find and alternative to convert the RATs to CSV or any format that allows me to do a join, after all.

I'm failing to identify a proper solution. Could anyone point me some clue? How can I convert rasters (kea format, in this case) to shape so that I can have the entire attribute table? Or can I export the RAT so that I can join it to a shapefile?

EDIT

Recently I made a progress. I'm able to store the column names (and eventually column values in python. But how to export/save them in a manageable format (cvs,txt,ascii) is still lacking solution.

>#!/usr/bin/env python
># Import python modules
>from rsgislib import rastergis
>from osgeo import gdal
># The output segments (clumps) image
>seg = "cur_segs_50_1000_UTM.kea"
>outfile = "test.txt"
>img= gdal.Open(seg)
>band=img.GetRasterBand(1)
>rat=band.GetDefaultRAT()
>icolcount=rat.GetColumnCount()
>cols=[]
>for icol in range(icolcount):
>cols.append(rat.GetNameOfCol(icol))


texts: Which projection can I use to see alignments on great circles?

texts:
I want to put on a QGIS map some historical sites to check if they are aligned.
If an alignement exists, the points should stay on a great circle (is this correct?).
I want to visually verify this using a map projection where orthodromic lines are represented as straight lines. In a second step I would draw two boundary lines and calculate bearings from a determinate point.

For example, here is a snapshot of a QGIS map with the points of "The St Michael-Apollo 'Axis'" as described by Lucien Richer, 1977.

Map boundaries: (-11E, 32N) (36E, 52N). EPSG 3785 projection (popular visualization / Mercator) .

My questions:

What would be the most appropriate map projection to get great circles as straight lines?
Do I need a custom projection? If yes, what are the parameters to use in QGIS?

Any suggestions would be helpful.


texts: How to connect Microstation V8i to Postgresql?

texts: I have been searching the web trying to find out if any one has connect Microstation to Postgresql, and to my surprise there has not been alot of chatter on ways to do this.. When looking at the database connection window in V8i I see no postgresql connection. Does anyone know a way to make this connection?


texts:
I'm curious if anyone has been able to change one individual value in a raster file?  I'm using Arc 10.2 and trying to change the value 99 to -2 for 144 files using model builder.  So far, I've only been successful using reclassify tool that requires me to manually eliminate the remaining value points.

Any suggestions?


texts: Getting the largest rectangle in the R-tree index of a PostGIS table?

texts:
Is it possible to use SQL or PL/pgsql to get the largest rectangle in the R-tree index (or other types of indices) of a table, so that it can be used as Bounding box for PostGIS table?

If possible, this would be much faster than going through all shapes and aggregating the bounding box (as suggested in the answers of Bounding box for PostGIS table).

(By comparison, shapefiles have bounding rectangles pre-computed. It would be nice if PostGIS can have this too in some form.)


texts: How To Use Tiled Map Services with Different Spatial References in ArcGIS API for JavaScript?

texts:
I have the following code fragment in an esri page to try and change basemaps, the following code works, unless the next basemap layer uses a different spatialreference.

if (baseMap)
  _map.removeLayer(baseMap);
baseMap = new esri.layers.ArcGISTiledMapServiceLayer(mapSource);
_map.addLayer(baseMap);

It seems like my only option is to destroy the map instance and start from scratch every time the user desires to change the basemap. Am I going about this the wrong way... anyone have any insight? I have seen samples of the basemapgallerywidget, but haven't used it.

So the basic question is: To change basemaps,not knowing the future basemap spatial reference, is my only option to destroy the map instance and re-create it?


texts: Modifying polygons to be more rectangular using PyQGIS

texts:
I have some buildings that were quickly digitized from an aerial image.  As you can see most of the polygons are not rectangular.
How can I modify the polygons to make them more rectangular?
Is there some equation that can solve a rectangle from a group of points?
I do not want to re-digitize the polygons.  The polygons are supposed to represent buildings.  I am using QGIS.
I tried code from @Matte but resulted in


texts:
I set up WMS server using MapServer, and the following requests work:

http://localhost:8081/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/example2_wms.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetCapabilities

http://localhost:8081/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/example2_wms.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=hydro&STYLES=&SRS=EPSG:4326&BBOX=-97.239,41.6198,-82.1229,49.3856&WIDTH=400&HEIGHT=300&FORMAT=image/png

I also tried to use the WMS that I serve in a map file, added the WMS in the following way:

LAYER
  NAME         hydro
  # PROCESSING   "LOAD_WHOLE_IMAGE=TRUE"
  TYPE         RASTER
  STATUS       ON
  CONNECTIONTYPE WMS
  CONNECTION "http://localhost:8081/cgi-bin/mapserv.exe?"

  METADATA
    "wms_srs" "EPSG:4326"
    "wms_name" "hydro"
    "wms_server_version" "1.1.1"
    "wms_format" "image/jpeg"
  END
END

Then requested it with the following URL: http://localhost:8081/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/example1-7.map&layer=hydro&mode=map

However, I only see a blank image.

Looks like this is the way, see http://mapserver.org/ogc/wms_client.html#wms-client & http://demo.mapserver.org/cgi-bin/mapserv?map=/osgeo/mapserver/tutorial/htdocs/example1-7.map&layers=states+states_label+modis_jpl&mode=map. Here modis_jpl is not visible too, but in the map file its status is off, this can be the reason.

Any idea on what I am doing wrong here?


texts:
I have three polygon shapefiles that each represent water (lakes, ponds, rivers) for three different time periods (1990, 2000, 2010). Each shapefile has over 1 million polygons/attributes and were generated from rasters, using a pretty long process in Arcpy (rough summary of the process: region grouping using all 8 neighbors, raster to polygon conversion, dissolve along the gridcode created by region grouping).

What I want to do is see how the water bodies have changed across these time periods, in particular, seeing how their areas have changed. In some cases, a lake/pond that existed in 1990 (for example) could have partially dried up and turned into 2 lake/ponds, or likewise, two lakes/ponds from 1990 could now be one polygon. The output doesn't have to be a shapefile, I just would (eventually) like a list of all the polygons represented in all of the shapefiles, and their corresponding areas for each time period.

So far, what I tried to do is a Spatial Join between the 1990 and 2000 shapefiles, then using that result (90_00_join.shp) in another Spatial Join with the 2010 shapefile (result is 90_00_10_join.shp). I wanted to capture the cases I mentioned above, so I used the Join one to many option. I also want all polygons from each time period to be present in the final output. So if a polygon (a.k.a. lakes/ponds) that existed in 1990 and/or 2000 but doesn't exist in 2010, I still want it in the final output. If a polygon didn't exist in earlier years but does in 2010, I also want it in the final output.

However, Spatial Join only allows for me to "Keep All Target Features." As a result, my final Spatial Join output has the same number of attributes as the 1990 input shapefile (the one that was used as the Target Feature in making 90_00_join.shp), although there are about 60,000 more attributes in the 2000 shapefile. So clearly these attributes from 2000 are getting left behind, I'm assuming because they are the 'Join feature' and not the Target feature. The same then happens after I join 90_00_join.shp with the 2010 shapefile.

Is there a way to get what I want using Spatial Join? If not, what about in python/arcpy by looping through the features and getting the area of each polygon and the polygons from the other shapefiles that intersect each feature? If not in Arc, what about QGIS or some other open source GIS software? The thing to consider is the number of polygons each shapefile has (about 1.2 - 1.3 million polygons).

EDIT: I should have mentioned that the problem I'm having is that the overlapping polygons do not have the same ID. So there is no way to compare the polygons across the years by just looking at the attribute tables. What I would like is to eventually have a table with ALL the polygons that existed throughout each time period, and their areas for each of the time period. This could normally be achieved using a spatial join, however what I can't figure out is how to get the output to include all polygons, not just those from the Target feature.


texts:
I want to to use an array to create 3 buffer files with different buffer distances. The code should produce 3 buffer files in one step. Additionally, the code below runs already after using an ArcGIS Add-in, but I have no idea on how to build in the array and to save this 3 different buffer files.

Can somebody please point me in right direction, as I have no idea where to start? My system is ArcGIS10.2, VB.net, Visual Studio2012

Here is the code without array:

Imports ESRI.ArcGIS.Geoprocessor
Imports ESRI.ArcGIS.Framework
Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI
Imports ESRI.ArcGIS.AnalysisTools

Public Class BufferButton1
    Inherits ESRI.ArcGIS.Desktop.AddIns.Button

  Public Sub New()
  End Sub

  Protected Overrides Sub OnClick()

    Dim m_application = My.ArcMap.Application
    Dim mxDocument As IMxDocument = m_application.Document
    Dim map As IMap = mxDocument.FocusMap

    Dim strInputFC As String
    Dim strOutputFC As String
    Dim intOutputBufferSize As Integer

    Dim gp As ESRI.ArcGIS.Geoprocessor.Geoprocessor = New ESRI.ArcGIS.Geoprocessor.Geoprocessor()
    gp.OverwriteOutput = 1

    strInputFC = map.Layer(0).Name
    strOutputFC = "E:\" ' in this path I would like to save the 3 buffer files
    intOutputBufferSize = 1000 ' here I would like to handle three buffer sizes eg 1000, 2000 and 3000

    Dim bufferTool As ESRI.ArcGIS.AnalysisTools.Buffer = New ESRI.ArcGIS.AnalysisTools.Buffer()
    bufferTool.in_features = strInputFC
    bufferTool.out_feature_class = strOutputFC
    bufferTool.buffer_distance_or_field = intOutputBufferSize
    gp.Execute(bufferTool, Nothing)
    My.ArcMap.Application.CurrentTool = Nothing
  End Sub

  Protected Overrides Sub OnUpdate()
      Enabled = My.ArcMap.Application IsNot Nothing
  End Sub
End Class


texts: Export shapefile script only works every other time, other times only exports XML

texts:
I am working on a simple python script that will check for the existence of a shapefile, remove that shapefile if it is found, then export an SDE feature class to a shapefile. The reason for doing it this way, is if the shapefile already exists, then when the tool exports the feature class to a shapefile it will add "_1" to the name, and so on.

I am having a problem with the tool when the script identifies that the shapefile exists in the user-defined folder. If the shapefile does not exist in the folder, the tool works flawlessly, however, if it does exists, I only end up with a .shp.XML instead of the actual shapefile.

Below is the script I am using. I have removed my actual database connection filepath and replaced it with the text "SDE database connection." This script prompts the user to choose an output folder as a parameter before initiating the script.

Why am I only getting the .shp.XML file extension when the script identifies that the shapefile already exists in the user-defined folder? Should I be using a different method to check for and delete the old shapefile?

import arcpy, os

# VARIABLES
Stormwater_Points = SDE database connection
Output_Location = arcpy.GetParameterAsText(0)

# REMOVE OLD SHAPEFILE
Old_Shapefile = os.path.join(Output_Location, 'Stormwater_Points.shp')
if os.path.exists(Old_Shapefile):
    os.remove(Old_Shapefile)
else:
    arcpy.AddMessage('Old shapefile not found in this directory, exporting new shapefile to selected directory')

# EXPORT NEW SHAPEFILE
arcpy.FeatureClassToShapefile_conversion(Stormwater_Points, Output_Location)
arcpy.AddMessage('Done')


texts: Networkx shortest paths for two point layers

texts:
I have a task to know the best way from houses to subway station.

There are three types of data on a picture below:

orange dots are the entrances to subway stations
violet dots are centerspoints of houses
blue lines is a road network (let it firstly contain all types of
road data)

The goal is to determine the shortest path from every house point to the closest subway entry and build a polyline layer.
I was thinking of using arcpy library for that but recently heard about networkx module, but right now I know only how to load data in graph, however, I cannot figure out how to look for desirable shortest path and what to do next in general.

Here is my draft for loading data to nx.Multigraph():

import arcpy
import networkx as nx
import numpy as np

line_fc = 'lines'
points_hs = 'house_points'

graph = nx.MultiGraph()
with arcpy.da.SearchCursor(line_fc, 'SHAPE@') as sc:
    for row in sc:
        line_geom = row[0].getPart(0)
        list_geom_cors = np.array([(n.X, n.Y) for n in line_geom])
        graph.add_edges_from(list_geom_cors)

with arcpy.da.SearchCursor(points_hs, 'SHAPE@XY') as sc:
    for row in sc:
        graph.add_node(row[0])

How to perform this task as correctly as possible?


texts: Short-term utility of ArcGIS/GIS certifications?

texts: (Correct me if this belongs in meta) For undergraduate students considering a career using GIS but who will not be graduating with a degree in it, is it a worthwhile investment of time to pursue a GIS certificate/certification? If so, what criteria are important to consider when selecting a test/program?


texts: Union, merge, intersection, join layers - logical difference and clear explanation for QGIS

texts:
I note that there are answers available on the difference between join and merge (e.g., Difference between 'join' and 'merge'?) that link to ArcGIS resources as a reference. However, I think that this problem is more complicated and I'm seeking a sleek summary applicable to QGIS.

I have a large data set that I integrated by running multiple unions between layers (e.g., "streams" U "amphibians" U "species at risk" U "ungulate winter range", etc...). The data set became too large to handle, so I had to break it into pieces and rebuild. This process has made me think more critically about the various union, merge, and join tools.

My goal is to identify what environmental features intersect with a study area that primarily consists of urban roads. I have four final layers that resulted from multiple rounds of unionizing layers of different environmental features:

P1P1
P1P2
P2P1
P2P2

Each has of the four layers has the same features within their respective data tables. I tried unionize these layers (e.g., P1P1 U P1P2 = P112, P2P1 U P2P2 = P221, and then P112 U P221 = Final) but I run into computation issues (see here).

As a work-around I thought I would run an intersection (∩) for each of the 4 layers and my study area (roads), which would reduce the computational power required. I thought I could then I merge the different resultant layers to get the result I'm after.

P1P1 ∩ "study area" = P1P1∩
P1P2 ∩ "study area" = P1P2∩...and so on...
Merge P1P1∩ with P1P2∩ = a, then
Merge P2P1∩ with P2P1∩ = b...and so on...and then
Merge a + b + c + d = Final result.

That was my plan. However, there is an issue with the data tables. For example, if I have a salamander and a frog in an overlapping area, then the merged data table will only report one or the other. It does not append the information to create the following:

The data table includes only the salamander OR the frog, but not both. Now, for comparison I ran merge, union (QGIS vector overlay), polygon union (SAGA), join attributes by location.

The merge creates a data table with 2079 rows and seems to be the correct output, but it only gives salamander OR frog and not (salamander, frog) where these spatially overlap.
The union (QGIS vector overlay) gives 2906 rows with separate fields for Amphibian and Amphibian_2, the polygon union (SAGA) gives 4979 rows with separate fields for Amphibian and Amphibian_1.
The join gives 2811 rows with separate fields for Amphibian and Amphibian_2.

I made this comparison because I thought there might be a way to append fields from the union data tables; while I can think of a way to do this with the field calculator, I would like to understand the logic prior to advancing in this direction. I was also surprised that polygon union (SAGA) did not give the same result as union (QGIS vector overlay), which makes me concerned and wonder about the steps that I followed to build the initial layers.

I have also tried to solve this problem by working with the rstats package sf using the geos_combine function. This creates another type of output that looks more like the merge function in the QGIS.

Is there a single illustrative resource that could help me to understand the difference in logic as it applies to the QGIS functions (union, merge, join, intersection) as this relates to the data tables as described in this example; same goes for r-stats sf package? Should I be seeking an Append function as exists in ArcGIS? I note that there is the plugin called "Append Features to Layer", but I want to clearly understand the logic to manage my data appropriately to bring this information together logically and with the correct final answer.


texts: Satellite/Aerial Imagery for Nepal for mapping in OpenStreetMap

texts: I was trying to map places in Nepal using satellite imagery but there is no high resolution imagery available at zoom level 15-16 . I was looking to map places like Illam, Dhankuta. Is there a free source of aerial imagery for these places so that i can use it with JOSM to map these places in openstreetmap?


texts: How can python-GDAL installed by Kyngesburye be recognized by pip?

texts:
I am using macOS 10.12 and Python 2.7 installed by homebrew.
As recommended, I successfully installed GDAL-1.11.4 by Kyngesburye's Framework.

When I want to install another package which dependents on "GDAL>1.9,<2.0" using "pip install *.whl" or "python setup.py install", the pip do not recognize the installed GDAL. And of course, the pip will download and install gdal automatically.

Then, error occurs with " File "", line 4, in
setuptools.sandbox.UnpickleableException: gdal_config_error(OSError(2, 'No such file or directory'),)"

Besides, when I edit my code using PyCharm, the GDAL is also not identified thought the code can run successfully.

I wonder how can I fix this problem?


texts: Memory management in offline routing with n waypoints?

texts:
I am trying to build an offline routing system with a map database on a portable Android device. Since I am still doing some research, I am stuck on what I think is a fundamental problem.

To find the best possible route between 2 points, one needs to load the road network of that area into memory and implement an algorithm which finds the shortest path.

If the two points are very far away from each other, the amount of road data that needs to be searched could be massive and impossible to be loaded in memory of a portable device running a system such as Android. The problem gets worse with even more waypoints.

Yet there are a lot of offline apps that perform such a thing. So it seems to me that my thought process is wrong.

Do these kinds of apps load the absolute necessary memory on the fly, discard and reload others for every road junction?

That seems such a slow process to me.


texts: Getting the 'memory:' layers after grass7:v.clean processing using PyQGIS?

texts:
To repair geometries errors, I'm now trying processing (grass7) using PyQGIS.
For one polygon layer (QgsVectorLayer: 'hn_pg'), I wrote the below script and run the grass7:v.clean algorithm. The processing completed without any problems, but the cleaned layer could not appear to the current active layers.

output = processing.runAndLoadResults(
            "grass7:v.clean",
           {'input': hn_pg,
            'type':[0,1,2,3,4,5,6],
            'tool':[1],   # snap
            'threshold':0.5,
            '-b':False,
            '-c':False,
            'GRASS_SNAP_TOLERANCE_PARAMETER':1,
            'GRASS_REGION_PARAMETER': "%f, %f, %f, %f"%hn_pg.extent().xMinimum(),hn_pg.extent().xMaximum(), hn_pg.extent().yMinimum(), hn_pg.extent().yMaximum()),
            'GRASS_MIN_AREA_PARAMETER':0.0001,
            'GRASS_OUTPUT_TYPE_PARAMETER':0,
            'output':'memory:',
            'GRASS_VECTOR_DSCO' : '',
            'GRASS_VECTOR_EXPORT_NOCAT' : False,
            'GRASS_VECTOR_LCO' : '',
            'error' : 'memory:'})

QgsProcessingUtils.tempFolder() function provided the path of the completed file(.gpkg) and I also found the outputted file in that file. So, I  tried QgsProject.instance().addMapLayer(output['output']) but resulted in no reaction.

How Can I get the cleaned object in Qgis active layers?


texts: Alter and update a shapefile with a field from an ASCII file using ogr

texts:
I have a vector file "map.shp" with a field ID.
I have an attribute table (data.txt) with the following fields: ID, value1, value2. The field ID is the same in both files.

I want to update map.shp by adding a new column and get the values from an attribute table ("data.txt").

How I could that with ogr in command line?

Since I don't have 2 shapefiles (vector layers).


texts: Donut shape related to concentric circles in QGIS

texts:
I would like to get an area between 10 miles and 20 miles from the center of counties. To do that, I wrote the code and got the following image. But, as you can see, there are some empty points, as the processing of "difference" clip the original 10 miles areas. I would like to get a donut shape. Any idea?
# Get centers of counties related to budgets
result = processing.run("native:meancoordinates", {
                        'INPUT': mem_layer,
                        'WEIGHT':None,
                        'UID':'fips_budgets',
                        'OUTPUT':'memory:'})
center = result['OUTPUT']

# Get buffers related to center from 10 miles
result = processing.run("native:buffer", {
                        'INPUT': center,
                        'DISTANCE':0.144985991,
                        'SEGMENTS':5,
                        'END_CAP_STYLE':0,
                        'JOIN_STYLE':0,
                        'MITER_LIMIT':2,
                        'DISSOLVE':False,
                        'OUTPUT':'memory:'})
buffer = result['OUTPUT']

# Get buffers related to center from 20 miles
result = processing.run("native:buffer", {
                        'INPUT': center,
                        'DISTANCE':0.289971982,
                        'SEGMENTS':5,
                        'END_CAP_STYLE':0,
                        'JOIN_STYLE':0,
                        'MITER_LIMIT':2,
                        'DISSOLVE':False,
                        'OUTPUT':'memory:'})
buffer2 = result['OUTPUT']

# Get 20-10 miles areas
result = processing.run("native:difference", {
                        'INPUT': buffer2,
                        'OVERLAY': buffer,
                        'OUTPUT':'memory:'})

buffer3 = result['OUTPUT']

QgsProject.instance().addMapLayer(buffer3)

UPDATE
I can get the result. Thanks! Actually, each center point has an ID (fips)[picture 2]. I would like to pass this information to the donut layer because I would like to define counties related to the center of the county. To do this, I am planning to use "intersection" such as the following code. Then, I would like to get information "This centered county (A) may affect the following counties (B, C, D) within the donut area". But, I could not get this information like [picture 3] because of the lack of the information.

Donuts = vl

# Get intersection of buffer and counties
result = processing.run("native:intersection", {
                        'INPUT': Donuts,
                        'OVERLAY': mem_layer,
                        'INPUT_FIELDS':[],
                        'OVERLAY_FIELDS':[],
                        'OVERLAY_FIELDS_PREFIX':'',
                        'OUTPUT':'memory:'})
intersection = result['OUTPUT']

QgsProject.instance().addMapLayer(intersection)


texts: I just recently started learning QGIS. I have a csv file with 10,000,000 entries detailing every yellow taxi trip taken in New York City in December of 2015. I want to create a trip distribution matrix with this data. I began by plotting the first 10,000 entries' pickup and drop-off longitude and latitude points for a total of 20,000 points. I then created a vector grid (.5 x .5 miles) in order to function as the zones for my trip distribution matrix. Then I tried the count points in polygon function in the vector analysis menu but the count comes up 0 in every row. Any idea what I am doing wrong?


texts:
I am trying to reconstruct some landscape features using the buffer option in QGIS. But, how can I modify it so it does an inward buffer, not outward?

Nothing really happens when I specify a negative value at the buffer distance. There are no buffers at all even though a new layer is shown.


texts: Searching many feature datasets for ArcSDE feature classes using field filters in ArcCatalog?

texts:
I have ~30 feature datasets and ~100 feature classes spread throughout the feature datasets.  I would like to search for a feature class by name among the feature datasets.  Also, I would like to query out of all the feature classes those records that have "Date_Mod" from someDate to someOtherDate.

This is all in ArcGIS SDE 10.0 on sql server 2008.


texts: Iterating raster and extract multi value of point files using ArcGIS Desktop?

texts:
I have several raster dataset (<8000 dataset) and one of point shapefile (contain 10 points).

I have created a model in ModelBuilder and extract it using multivalues to points feature on ArcGIS.

Can you help me to correct my model?

Are there any possibilities to create a field name on shapefile with same name as raster file?

For example I have raster file "201702022300.tiff" and on "buoy_sawu.shp" I will have an output fieldname (using multivalue to point tools) a "201702022300" field name?


texts: Negotiating ArcPy schema locks in script

texts:
How do I avoid the 000464 (Cannot get exclusive schema lock) error in this code caused by the last line?
import arcpy
import os

geodatabase = 'C:/Data/testgdb.gdb'
fc_name = 'test_feature_class'

fields = [
    ['name', 'TEXT', 'Name', 32],
    ['other', 'TEXT', 'Other', 32]
]

if not arcpy.Exists(geodatabase):
    gdb_path, gdb_filename = os.path.split(geodatabase)
    print(gdb_path, gdb_filename)
    gdb_name, gdb_extension = os.path.splitext(gdb_filename)
    print(gdb_name, gdb_extension)
    gdb = arcpy.CreateFileGDB_management(gdb_path, gdb_name)

if not arcpy.Exists(os.path.join(geodatabase, fc_name)):
    fc = arcpy.CreateFeatureclass_management(geodatabase, fc_name, 'POINT', spatial_reference=arcpy.SpatialReference(4326))
    ## THE BELOW LINE TRIGGERS ERROR 000464: Cannot get exclusive schema lock.  Either being edited or in use by another application or service. ##
    arcpy.AddFields_management(fc, fields)

Interestingly, the entire code executes without issue if I delete the entire geodatabase between runs (not an option for production, just an observation).
If I remove the offending arcpy.AddFields_management() method call, the code will create the feature class in the existing geodatabase without error (and without schema).  But I cannot call AddFields_management() on the feature class created in the existing database without having it raise a schema lock error.


texts: How to convert micostation v8 dgn to IFS in FME?

texts:
I have a 3D model saved in dgn format as V8. This may be initially created in Bentley AECOsim. I can't able to convert it to IFC format?

I added reader as dgn(v8) and writer as IFC in FME workbench. It converts some feature in a particular layer but not all. In log it seems it successfully read all feature but couldn't write all. What ever the feature got write is not view-able in FMEinspector.

Can anyone point in good direction?


texts: How can I have ArcGIS not automatically add raster clip output from Python window to the current map

texts:
I have a collection of rasters that I'm writing code to automate clipping based on a user defined shapefile. I am currently doing this in the Python window using the following as the main processing function:

arcpy.Clip_management( inputRasterFilename, "#", outputRasterFilename, "shapefile" )

Basically I want to read the raster, clip it, and have it saved to outputRasterFilename for later use.

This is exactly what happens except that ArcGIS also adds the result of clip operation to the current map that I have open. Is there a way to disable this?


texts: Where is "Nearest neighbour analysis" tool in pyqgis documentation?

texts:
I have a group of point features for which I want to run the Nearest neighbour analysis using python. The problem is that I cannot find that function in the documentation.

The closest I got was:
https://qgis.org/pyqgis/3.0/core/other/QgsSpatialIndex.html?highlight=nearest#qgis.core.QgsSpatialIndex.nearestNeighbor

However, this seems to work on a point level? I want to run the analysis per feature like you would do using Vector Analysis->Nearest neighbour analysis.

Could someone pin me this in the documentation, or suggest an alternative?


texts: Replace multiple attribute field headers in PyQGIS

texts:
I am new to pyqgis plugins with QT, I am writing code to replace multiple headers in the attribute field of a layer. I keep getting an error saying

  'int' object has no attribute 'pendingFields'

What can I do?

Here is my code. The code accepts the current shapefile and renames the headers:

def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    #This activates the layer on the table of content
    self.dlg.layercombo.clear()

    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layerType = layer.type()
        if layerType == QgsMapLayer.VectorLayer:
            layer_list.append(layer.name())
    self.dlg.layercombo.addItems(layer_list)

    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        self.selectedlayerIndex = self.dlg.layercombo.currentIndex()
        for field in self.selectedlayerIndex.pendingFields():
            if field.name() == 'CC':
                with edit(self.selectedlayerIndex):
                    idx = self.selectedlayerIndex.fieldNameIndex(field.name())
                    self.selectedlayerIndex.renameAttribute(idx, 'DD')

            if field.name() == 'AA':
                with edit(self.selectedlayerIndex):
                    idx = self.selectedlayerIndex.fieldNameIndex(field.name())
                    self.selectedlayerIndex.renameAttribute(idx, 'BB')


texts: Checking if a geocoordinate point is land or ocean?

texts:
I got the corner geo coordinates from a xml product file of a Sentinel-1 GRD SAR image. When looping through the pixels of the image I wish to check whether it represents land or sea.

I figured I could calculate the corresponding longtitude and latitude by using the corner longtitudes and latitudes and do the check. Looking into basemap from mpl_toolkits, i found an attribute .is_land(x, y) but this require a map projection which slows it down, which I want to avoid.

Is there an efficient and fast way to check whether a longtitude and latitude is land or ocean in python 3.x?


texts:
I create wfs vector layers in loop:

 for(var i=0; i<vectors.length;i++)  {
                 var styles = {
                            'Point': new ol.style.Style({
                                image: new ol.style.Circle({
                                    radius: vectors[i].pointRadius,
                                    stroke: new ol.style.Stroke({
                                        color: vectors[i].pointStrokeColor
                                    }),
                                    fill: new ol.style.Fill({
                                        color: vectors[i].pointFillColor
                                    })
                                })
                            }),
                            'MultiPoint': new ol.style.Style({
                                image: new ol.style.Circle({
                                    radius: vectors[i].pointRadius,
                                    stroke: new ol.style.Stroke({
                                        color: vectors[i].pointStrokeColor
                                    }),
                                    fill: new ol.style.Fill({
                                        color: vectors[i].pointFillColor
                                    })
                                })
                            }),
                            'LineString': new ol.style.Style({
                                stroke: new ol.style.Stroke({
                                    color: vectors[i].lineColor,
                                    width: vectors[i].lineWidth
                                })
                            }),
                            'MultiLineString': new ol.style.Style({
                                stroke: new ol.style.Stroke({
                                    color: vectors[i].multiLineColor,
                                    width: vectors[i].multiLineWidth
                                })
                            }),
                            'Polygon': new ol.style.Style({
                                stroke: new ol.style.Stroke({
                                    color: vectors[i].polygonStrokeColor,
                                    lineDash: [vectors[i].polygonStrokeLineDash],
                                    width: vectors[i].polygonStrokeWidth
                                }),
                                fill: new ol.style.Fill({
                                    color:vectors[i].polygonFillColor
                                })
                            }),
                            'MultiPolygon': new ol.style.Style({
                                stroke: new ol.style.Stroke({
                                    color: vectors[i].multiPolygonStrokeColor,
                                    width: vectors[i].multiPolygonStrokeWidth,
                                }),
                                fill: new ol.style.Fill({
                                    color: vectors[i].multiPolygonFillColor
                                })
                            })
                        };
                        var styleFunction = function (feature) {
                            return styles[feature.getGeometry().getType()];
                        }
                        overlays.push(
                              new ol.layer.Vector({
                                    visible:vectors[i].isDefault,
                                    title: vectors[i].name,
                                    style: styleFunction,
                                    source:new ol.source.Vector({
                                       format:new ol.format.GeoJSON({
                                         featureProjection: "EPSG:3857"}),
                                    url:vectors[i].url
                            })}));
}

After that I add overlays array of vector to map and styles don't work. All vectors have style from last vector in array. How make difference vector style?


texts: cannot append a featureCollection in typescript

texts:
I have the following interfaces:

export interface Feature<T extends GeometryObject> extends GeoJsonObject
{
    type: 'Feature'
    geometry: T;
    properties: any;
    id?: string;
}

export interface FeatureCollection<T extends GeometryObject> extends GeoJsonObject
{
    type: 'FeatureCollection'
    features: Feature<T>[];
}

The I tried to append the following objects:

let geoJsonFeatureCollection:FeatureCollection<any> = { type: 'FeatureCollection',

features: []
};

console.log(geoJsonFeatureCollection);

for (var key in data){
    if (data[key].Volume>0){
        let Feature:Feature<any>;
        Feature['type'] = 'Feature';
        Feature['geometry'] = {'type': 'Point',
        'coordinates': [data[key].YO, data[key].XO],
    };
    Feature['properties'] =  {'origin_city': data[key].Origin,
    'origin_lat' : data[key].YO,
    'origin_lon' : data[key].XO,
    'destination_city' : data[key].Destination,
    'destination_lat' : data[key].Y,
    'destination_lon' : data[key].X,
    'volume' : data[key].Volume

};
    geoJsonFeatureCollection['Features'].push(Feature);
}
}

but looks like its not working ..


texts: Geoserver ImageMosaic BackgroundValue Ignored for non-RGB GeoTIFF

texts:
I have a GeoServer (2.10.0) instance into which I have loaded some GeoTiff files, one per timestep, using the ImageMosaic plugin.  I can query them with standard WMS requests with the parameter TIME=yyyy-mm-dd.  The relevant data is quite sparse and I have created an SLD style which renders the background data as transparent.

What I would like is to pass a string such as TIME=2013-01-01/2015-01-01 and have all of the images which fall within that time range returned and layered on top of one another.  What happens is that I get the latest image in that time range returned.  Is what I want to do possible with GeoServer, and if so how can I enable it?

EDIT: After further investigation it appears the GeoServer is layering the data, but:

The background data value in the GeoTIFF is zero and NullValues in the Coverage Band details is set to zero, as is BackgroundValues
Despite this, large regions of background data from the latest time are obscuring all the data from the earlier times
The regions of no data, once layered and obscuring older data are then correctly rendered as transparent (as defined in the SLD style)
I can't set the InputTransparentColor (or OutputTransparentColor), I think because the GeoTIFF is not RGB.  If I try, I get:

org.geotools.data.DataSourceException: Unable to create this mosaic
java.lang.IllegalStateException: Unsupported data type.

So I think the problem is now that when using the ImageMosaic plugin, I cannot configure the value 0 to correspond to missing data properly.

I have tested all 16 combinations of:

MergeBehaviour: FLAT or STACK
BackgroundValues: 0 or empty
Coverage Band Null Values: 0 or empty
Coverage Band minRange: 0 or 1 (the values in the GeoTiff range from 0-4, with 0 representing no data)

None of these 16 combinations work.  Any ideas?


texts: What are the white grids at 9.6953038715817,98.576672019046 in Google Maps?

texts:
A friend asked me to help interpret the white grids at 9.6953038715817,98.576672019046 in Google Maps. In satellite image the grids are connected to rivers flowing into the Andaman sea. Most of the rivers are also covered with snow-white material. I guess the grids are for polluted water treatment. But these grids seem to lack patterns of water treatment plants.

Can anyone help interpret these grids? Clues are also appreciated!


texts: Define Terrain/Elevation textfile data format

texts:
I have a program that can visualize terrain/elevation data but will only recognize files in the below format. It has the extension .map but it is not recognized by QGIS or GDAL so I believe it may be another form of data.

When opened in a text file, it appears like this:

Example Map
XScale: 8.455
YScale: 10.447
XOrigin: 606319.184
YOrigin: 4056521.23
NumberOfCols: 91
NumberOfRows: 73
RangeUpper: 68.81
RangeLower: 30
0057.70 0058.18 0058.73 0059.29 0059.78 0060.18 0060.51 0060.81 0061.14 0061.46 0061.71 0061.81 0061.73 0061.52 0061.26 0061.07 0060.98 0060.98 0061.03 "...

The strings of data go on and on. Can anyone identify this format?

Would like to be able to convert to this from other elevation data types: GeoTIFF, DEM.


texts: How to get a text character based on a number position in the alphabet in a QGIS expression

texts:
This seems very simple but I can't find a function that will take a number between 1 and 26 and return the right letter.
1=A
2=B
3=C
...
Z=26


texts: Extract height from DTM into attribute table of point layer

texts:
I am using ArcGIS Desktop 10.4

I have an empty point layer (100 random points) and a dtm.

I was wondering how I would go about to get the height from the dtm underneath each point into the points attribute table?


texts:
I am trying to create some visualisations using bus journey data. I am looking to plot journeys along one long road, showing passenge journeys. The principle is the same as a radial flow map for flight paths for example. However the geographical area is much smaller as it only covers one road.

This is where I am running into problems. I have the data in StartX,Y and End X,Y format and then with each possible journey I have a total number.

I have tried plotting this as XY to Line using the various different Line Types (e.g. Geodesic). The result is a series of straight lines on top one another. Ideally I would like to produce a series of curved lines so that you can see the different journeys a bit more clearly.

I have tried to research this problem quite extensively online but the majority of cases are plotting data over a much wider geographical area.

Is there a way to produce a curve / arc between a start and end point?

I do not want to do this manually as I have large number of lines to plot.

I would like to use ArcMap.


texts: Exporting map containing raster layer with vector field symbology ends up with partially rendered symbology from ArcMap?

texts:
After exporting a map in ArcMap (ArcGIS 10.3 for Desktop) which contains several raster and vector layers including one raster layer resembling flow direction of ice, the output file (tried JPG, PNG and PDF so far) always contains lines in which the Vector Field Symbology is only partially rendered.

In this very example it seems like a pattern of four perfectly fine rendered symbols, followed by four lines which a partially rendered (repeating pattern).

Screenshot in ArcMap Layout View:

Crop of exported map:

Any ideas how to fix this?


texts:
I'm wondering if there is any suitable algorithm that can I use to calculate line complexity index (or indicator)? By "line complexity" I mean how many bends (direction changes) that line contains and how sharp they are.

I want to use that index to correct driving speed of road network segments. Lots of road segments in my dataset have speed defined according to road type, which is not even close to the real achievable speed. For example, mountain pass road has a speed limit of 90kph, but in real life, 60kph is more realistic due to curvature of that road.

I'm primarily interested in a solution for PostGIS/PostgreSQL.

The best solution i have come up so far is something like this:

st_area(st_polygonize(st_union(st_simplify(geom, 10), geom)) / st_length(geom)

It calculates how much space (area) per unit of length is between the original and simplified version of geometry. Bigger index means more complex shape, thus bigger speed penalty. Simplification threshold (10 in my example) was chosen by trial and error, also how it will be applied to adjust speed is not clear at this point.

Is this feasible, or is there a better approach?