texts: Proj4Leaflet with EPSG:2263 does not seem to project correctly

texts:
For some reason EPSG:2263 projections seem to be incorrect.  The tile server works properly at all zoom levels as configured.  However, the tiles are rendering at the incorrect latLong.

Origin was obtained from server's bounds: Only Top left (tl) works.
the transformation solution commented in code has no effect

Current Proj4Leaflet config:

var empireStateBldg = new L.latLng(40.7484444, -73.9878441),
// Bounds from server: Top left (tl) works as origin.
tl = L.point(700000.0,  440000.0),
br = L.point(1366666.6683464567, -4444.4455643044785),
tr = L.point(1366666.6683464567, 440000.0),
bl = L.point(700000.0, -4444.4455643044785);

var nycCrs = new L.Proj.CRS("EPSG:2263",
    "+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 " +
    "+lon_0=-74 +x_0=300000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 " +
    "+units=us-ft +no_defs",
    {
    resolutions: [
        434.0277788713911, 303.8194452099737, 222.22222278215222, 111.11111139107611, 55.555555695538054,
        27.777777847769027, 13.888888923884513, 6.944444461942257, 3.4722222309711284, 1.7361111154855642,
        0.8680555577427821, 0.43402777887139105, 0.21701388943569552, 0.10850694471784776
    ],
    origin: [tl.x, tl.y],

    //https://github.com/kartena/Proj4Leaflet/issues/4
    // (1, -x_origin, -1, y_origin)
    // transformation seems to have no effect
    transformation:new L.Transformation(1, -700000.0, -1, 440000.0),
    })

var map = L.map('map', {
    crs:nycCrs,
});

L.tileLayer.wms('http://maps.nyc.gov/geowebcache/service/wms/', {
    layers: 'dtm',
    format: 'image/png',
    maxZoom: nycCrs.options.resolutions.length,
    minZoom: 0,
    continuousWorld: true,
    tileSize: 512,
    worldCopyJump: false
}).addTo(map);

map.setView(empireStateBldg, 0);

L.grid().addTo(map);

L.marker(nycCrs.projection.unproject(tl)).addTo(map);
L.marker(nycCrs.projection.unproject(tr)).addTo(map);
L.marker(nycCrs.projection.unproject(bl)).addTo(map);
L.marker(nycCrs.projection.unproject(br)).addTo(map);
L.marker(empireStateBldg).addTo(map);

var popup = L.popup();

function onMapClick(e) {
    popup
        .setLatLng(e.latlng)
        .setContent("You clicked the map at " + e.latlng.toString())
        .openOn(map);
}

map.on('click', onMapClick);


texts: Editor Tool, Areas instead of Segments in ArcMap?

texts:
I'm studying this tutorial:

https://learn.arcgis.com/en/projects/get-started-with-arcmap/lessons/predict-the-impact-of-the-proposed-road.htm

But I can't complete the "Digitize the proposed road" section.

When I use the "Trace" and the "Straight Segment" tools, at the end of the drawing I have those "areas" and not just the segments that I would need (you can see the difference looking at picture at point 16 in the previous link):

Is there some option I should be turning off?


texts: Vector Tiles: do the zoom levels depend on the projection?

texts: For Mercator Projection I found a table of zoom levels and their corresponding scale. I suppose that the number of zoom levels and the corresponding scale change when using another projection for Map Tiles Services?


texts: Iso-Area jumping over parts without streets (QGIS QNEAT3 Plugin)

texts:
I am using the QGIS Plugin 'QNEAT3' (version 1.0.4) to calculate Iso-Areas. They are only relevant for me in one country, so I only have this country's street network. But the Iso-Areas are jumping over the parts without streets outside the country and then continue in another part of the country.
To show what I mean I provided a map with my result.
How can I prevent this and make the Iso-Areas stop at the border of the country?


texts: Selecting polygon if their centroid falls in another polygon using PostGIS?

texts:
I'm trying to select polygons which centroids fall in another polygon:

SELECT smaller_polygons.id
FROM big_polygon
LEFT JOIN big_polygon smaller_polygons
ON ST_Within(ST_Centroid(smaller_polygons.wkb_geometry),big_polygon.wkb_geometry)

I expect 33,44 to be selected but I get id of the big polygon.

What am I doing wrong?


texts: ReduceRegion(s) Google Earth Engine Python API

texts:
I'm trying to find the area of a classified land cover across the study area using the Google Earth Engine's Python API. I am trying to do this by mapping a .reduceRegion function across the image collection using the approach below.
# Salt marsh extents are images where pixels = 1 if salt marsh has been detected

salt_marsh_area = salt_marsh_extents.map(lambda image:
                                      image.reduceRegion(
                                          reducer=ee.Reducer.sum(),
                                          geometry=whole_estuary.geometry(),
                                          scale=30
                                      ))

However, when I try and export this to my Google Drive or use a .getInfo() to inspect the result I get an error of ee.ee_exception.EEException: Collection.map: A mapped algorithm must return a Feature or Image.
How do you convert this into a table format for export?


texts:
I have a  .csv file with several attribute fields. I would like to join this to a raster (.tif) file, using an 'ID' column from the CSV file, so that the raster gets all the attributes from the CSV file.

Here's the R code for reading TIFF and CSV:

library(raster)
library(tidyverse)
     # Read tif file with values that are long integer type
        tif_arc <- raster("tif_arc.tif")

        #Read csv file that has an attribute 'ID' the value of which = the raster value on the tif file
        csv_whole <- read_delim("csv_qgis.csv",delim = ",")

    glimpse(csv_whole)

Here's what the CSV
data looks like on running "glimpse()" command:

$ FOLIAGE_BIOMASS_PER_HA        <dbl> NA, NA, NA, NA, NA,...

$ BARK_BIOMASS_PER_HA           <dbl> NA, NA, NA, NA, NA,...

$ ID                            <dbl> 13453064, 11242315, 11150335,...

I know this join can be done in ArcGIS using 'ADD JOIN' tool, using attribute field 'ID' but can this be done in R?


texts: Free hosting of shapefiles to make an interactive map?

texts: What is a free site where I can upload shapefiles and make my own interactive map to share with others by a link?


texts:
How to return raster values based on polygon (shapefile) boundaries in QGIS?

Is this an elementary question, because I am finding it hard.......
Ideally I would like a neat workflow as I am likely to repeat the task, even if it includes using GRASS.

I want to take a remote sensed raster ( Landsat ), and load it up with a shapefile containing irregular polygons and return the avg, max min of the raster values for each polygon.

Ideally I would also like to clip out the polygons as separate files for each polygon, perhaps that is another post, but probably it is a related process, which I haven't been able to tackle yet.

And is QGIS the right platform to be tackling this in? I have loaded up the RASTERCALC Plugin which does give better functionality than the default raster calculator, so have made a little progress there.

Any thoughts appreciated


texts: Failure to Iteratively Calculate Zonal Statistics as Table ArcGIS

texts:
I want to calculate mean luminosity for a district using several rasters and a shapefile with administrative borders. Since I have many files, I want to do the operation iteratively using a model builder in ArcGIS. This how the model looks like:

My logic was as follows:

Take every raster from the folder `Lights_2013_2016'
Calculate zonal statistics
Save statistics as a .dbf file using the original name of the file

Yet, when I perform it, I get a mistake:

  Executing (Zonal Statistics as Table): ZonalStatisticsAsTable C:\Users\1\Desktop\Projects\East_Ukraine\Shapefiles\UKR_adm\UKR_adm2.shp NAME_2 D:\Lights_2013_2016\20130101.tif C:\Users\1\Documents\ArcGIS\Default.gdb\20130101.tif.dbf DATA ALL
  Start Time: Sun Nov 12 21:12:21 2017
  ERROR 999999: Error executing function.
  Failed to execute (Zonal Statistics as Table).
  Failed at Sun Nov 12 21:12:23 2017 (Elapsed Time: 1,43 seconds)

The operation works manually. The operation does not perform in the Model Builder.
Could you tell me, what did I do wrong and what should I do?


texts:
I want to split polygons with multilines, and merge results with the neighboring polygon.

I tried with this code :

--Polygon splitting
CREATE TABLE decoupe AS
SELECT rect.affectatio, (
    ST_Dump(
        ST_Split(rect.geom, st_union(line.geom)
                )
    )
).geom as geom
FROM total_avant_decoupe as rect,
ensemble as line
GROUP BY rect.geom, rect.affectatio;

--Merging with neighboring polygon
CREATE TABLE fusionvoisinmajoritaire AS
WITH tempo AS (
    SELECT DISTINCT ON (id) id, ST_Union(geom) as geom, affectatio
        FROM (
        SELECT DISTINCT ON (b.id) a.id, b.id AS idb, ST_Union(a.geom, b.geom) as geom, a.affectatio
        FROM (
            SELECT *
            FROM fusiondifference
            WHERE affectatio != 'non intersecte'
        ) a JOIN (
            SELECT *
            FROM fusiondifference
            WHERE affectatio = 'non intersecte'
        ) b ON ST_Intersects(a.geom, b.geom)
        WHERE ST_Length(
                ST_Intersection(a.geom, b.geom)
            ) > 0.05*ST_Perimeter(b.geom)
    ) AS sub
    GROUP BY id, affectatio
)
SELECT *
FROM tempo
UNION ALL
SELECT *
FROM (
    SELECT *
    FROM fusiondifference
    WHERE affectatio != 'non intersecte'
) AS sub
WHERE id NOT IN (
    SELECT id
    FROM tempo
);

But with my splitting code, I get holes between my polygons. So in many case, it can't be merge with my merging code...

I tried also by snapping with ST_SnapToGrid and ST_Snap with the neighboring polygon but it doesn't work, I keep holes.

How to split my polygons without create holes between them? I'm using PostGIS 2.5.3.

EDIT (13/12/2019) :

This is a schema of that I want to do :

I have 3 options :

The first is simple, just one adjacent polygon, no split needed
The second just needed to split polygon by single lines
The third is little more complex. I create lines with ST_ApproximateMedialAxis amongst others. This requires two steps for the external polygons who haven't neighboor.

Make lines (in red) is not a problem. But the holes creates by splitting...

EDIT 2 (16/12/2019):

For illustrate, when I merge all polygons in one, I obtain this type of polygon, that show the splitting code make holes :


texts: osm2pgrouting stuck (showing classes table created) when importing an osm file

texts:
I am trying to import an osm file using osm2pgrouting. while following the pgrouting workshop. (file is 700mb~)

im running ubuntu 12.04 on m1.large on amazon.

the same file was successfully upload to postgres with osm2pgsql.

i am using the following command:

osm2pgrouting -file "/home/ubuntu/israelmap/israel_and_palestine.osm"  -conf "/usr/share/osm2pgrouting/mapconfig.xml" -dbname routing -user postgres -clean

the proccess starts ok, but the way table will never have any data in it:

We need a way of type highway and class service
We need a way of type highway and class service
We DON'T need a way of type  and class
We DON'T need a way of type  and class
We DON'T need a way of type  and class
Split ways
NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index "nodes_pkey" for table

"nodes"
    Nodes table created
    NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index "ways_pkey" for table "ways"
    Ways table created
    Types table created
    Classes table created
create topology
NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 15 at EXECUTE statement
#########################
size of streets: 159747
size of splitted ways : 365597
finished

after a very long while (2 hrs) it sais finished but way table is empty!
how do i fix this?

i was able to load a routing table using osm2po. still no success with osm2pgrouting

Thanks for your help!


texts: After an update, reset the values with QGIS

texts:
I'm working on a project which is using QGIS/PgAdminIII and I would like to know how to reset values after updating a specific one.
I created a SQL function with PgAdminIII console to reset everything after an update but you need to 1/ edit the value, 2/ save, 3/ reopen the form to get all the values reseted, but I'd like to save time and avoid people to open/close the form.
In details, every time I change the value of "brood" I would like all the others to return empty, to be able to start again entering data and saving time without opening/closing windows' form.

I almost find a way by using Qt designer and configuring the connection between "brood" an the others items but it failed. But I could have done a mistake setting up the connections. Any clues ?
An other option I have is to create a python code and loading it in the properties of the layer but my abilities are too basics to create such a thing. Even though I searched in the web I didn't find any appropriated case but I found this command could help me, right ?
layer.changeAttributeValue(fid, fieldIndex, value)

will it be possible to find a way to reset these values after updating brood?


texts: Arcpy SpatialJoin performed in_memory not returning any rows

texts:
Background: I am writing a tool to automatically project a layer (vector or raster) in ArcGIS Pro to its appropriate UTM zone.
Approach: The entire UTM zone grid is stored in a geojson string inside the script.  From the string, create a temporary polygon feature class in the in_memory space, pulling only the EPSG code in as a field value.  Get the centroid of the raster or vector layer.  Run a spatial join between the centroid and the UTM grid.  Obtain the UTM zone EPSG code.  Run a Project() (if vector layer); run a ProjectRaster() (if raster layer); output the results to the project's default geodatabase.
Problem: arcpy.analysis.SpatialJoin() is not returning any rows, even though it looks like I've got correctly constructed feature classes in in_memory.
Code:
import arcpy
from arcpy.sa import *

# methods
def make_temp_utm_fc(gj_string):
    sr = arcpy.SpatialReference(4326)
    temp_fc = arcpy.management.CreateFeatureclass('in_memory', 'utmzones', 'POLYGON', spatial_reference=sr)
    fields = ['EPSG']  # list in case i want to add more attributes later
    for field in fields:
        arcpy.management.AddField(temp_fc, field, 'TEXT')

    with arcpy.da.InsertCursor(temp_fc, ['SHAPE@'] + fields) as irows:
        for feature in gj_string['features']:
            geom = arcpy.AsShape(feature['geometry'], False)
            irows.insertRow([geom] + [feature['properties']['EPSG']])

    return temp_fc

def get_extent_centroid(layer):
    desc = arcpy.Describe(layer)
    sr = desc.spatialReference
    if desc.dataType in ['FeatureLayer', 'RasterLayer']:
        extent = desc.extent
        array = arcpy.Array()
        array.add(extent.lowerLeft)
        array.add(extent.lowerRight)
        array.add(extent.upperRight)
        array.add(extent.upperLeft)
        array.add(extent.lowerLeft)
        pg = arcpy.Polygon(array, sr)

        return pg.centroid

# env
p = arcpy.mp.ArcGISProject('CURRENT')
arcpy.env.Workspace = p.defaultGeodatabase
arcpy.env.overwriteOutput = True

# utm zones (abbreviated for here)
utmgrid = {
    "type": "FeatureCollection",
    "name": "utmgrid",
    "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},
    "features": [
        {"type": "Feature", "properties": {"fid": 1.0, "SWLON": "-180", "SWLAT": "-80", "HEMISPHERE": "s", "ZONE": "01", "Zone_Hemi": "01,s", "EPSG": "32501"}, "geometry": {"type": "Polygon", "coordinates": [[[#the coord ring]]]}}]}

# do the work
inlyr = arcpy.GetParameter(0)
inlyr_desc = arcpy.Describe(inlyr)

utm_fc = make_temp_utm_fc(utmgrid)

center = get_extent_centroid(inlyr)
center_row = [(center.Y, center.X)]
center_fc = arcpy.management.CreateFeatureclass("in_memory", "layercenter", "POINT", spatial_reference = inlyr_desc.spatialReference)
with arcpy.da.InsertCursor(center_fc, ['SHAPE@XY']) as cursor:
    cursor.insertRow(center_row)

# this is the part that returns nothing
sj = arcpy.analysis.SpatialJoin(center_fc, utm_fc, "in_memory/sj")

with arcpy.da.SearchCursor("in_memory/sj", "EPSG") as cursor:
    for row in cursor:
        arcpy.AddMessage(row)  # returns None

Once I have the EPSG code from the containing UTM zone polygon, it's easy to do the reprojection so I haven't included that code.  Where am I going wrong doing the spatial join?
Currently I'm working exclusively with layers that carry a 4326 spatial reference, so both the centroid and UTM grid are in GCS_WGS_1984.  I've added AddMessage() calls throughout for debugging and it seems like it's successfully created the UTM grid and centroid as Polygon and Point feature classes respectively.


texts: Available snapping points for advanced digitizing?

texts:
In QGIS, the Advanced Digitizing tools have some snapping features that I don't see mentioned in the documentation. We know that we can digitize new segments that are perpendicular or parallel to existing segments, or we can lock the new segment to any angle. However, having done any of these, we can also set the length of the new segment by snapping to a variety of "reference points" on other existing segments and vertices. For example, we can snap:

To the intersection of an existing segment, OR an imaginary extension parallel to that segment (a blue reference line is shown in QGIS):
To a vertex on an existing segment that would be at right angles to our new segment (here, I've highlighted the reference vertex with a red circle, and added a dotted green line to show the imaginary perpendicular reference line):

My question is, can somebody fully explain how this feature works, or is it documented somewhere? Are there any other ways to snap to a reference vertex or segment, besides the examples I've given?


texts: Turning off visibility of layer using ArcObjects Java?

texts:
I'm making a geoprocessing tool with ArcObjects Java (IBasicGeoprocessingTool) and I need to make a layer not visible. I managed to get this layer as IFeatureLayer and I get an access to the map or activeView.
I tried that:

featureLayer.setVisible(false);
((IActiveView)map).refresh();

After that, the layer is not visible but it remains checked in the table of contents. Does anyone know how to access to the TOC or MapControl to uncheck it please?


texts: Calculate where a geodesic line crosses a meridian

texts:
Given a geodesic line between two points with known coordinates on a sphere, is there an easy way to calculate or approximate where this line crosses a specific meridian?

I want to calculate where the shortest path between P1 and P2 crosses a meridian. So far all the documentation I found was either above my math level or would take me ages to work through, GEODESICS ON AN ELLIPSOID - BESSEL'S METHOD for example.


texts: Overpass install fails with a broken pipe error

texts:
I've installed a Ubuntu box which will be my one of many overpass api servers.
Basically the installation is good, and I've imported many testing areas including andorra-latest.osm.bz2 and taiwan-latest.osm.bz2. I can't however import china-latest.osm.bz2.

Here is the script I use to import

//this one failed
../bin/init_osm3s.sh china-latest.osm.bz2 "db/" "../" --meta=no

//this executes with no problem at all
../bin/init_osm3s.sh andorra-latest.osm.bz2 "db/" "../" --meta=no

It takes about 10 seconds to import andorra file with no problem at all. all db files created. But when I tried it with china file, I got this error. I don't have any clue at this moment.

../bin/init_osm3s.sh china-latest.osm.bz2  "db/" "../" --meta=no
Reading XML file ... elapsed node 2819628408. ../bin/init_osm3s.sh: line 44:
12709 Broken pipe             bunzip2 < $PLANET_FILE
 12710 Killed                  | $EXEC_DIR/bin/update_database --db-
dir=$DB_DIR/ $META $COMPRESSION

I downloaded those files using the following links

http://download.geofabrik.de/asia/china-latest.osm.bz2

http://download.geofabrik.de/europe/andorra-latest.osm.bz2

Need your help, please, if you know why or at least theories of how this happened.

P.S. init_osm3s.sh script

if [[ -z $3 ]]; then
{
  echo "Usage:  $0  Planet_File  Database_Dir  Executable_Dir  [--meta]"
  echo "        where"
  echo "    Planet_File is the filename and path of the compressed planet     file, including .bz2,"
  echo "    Database_Dir is the directory the database should go into, and"
  echo "    Executable_Dir is the directory that contains the executable     update_database."
  echo "    Add --meta in the end if you want to use meta data."
  exit 0
}; fi

PLANET_FILE="$1"
DB_DIR="$2"
EXEC_DIR="$3"
META="$4"
COMPRESSION="$5"

if [[ ! -s $PLANET_FILE ]]; then
{
  echo "File $PLANET_FILE doesn't exist."
  exit 1
}; fi
mkdir -p "$DB_DIR/"
bunzip2 <$PLANET_FILE | $EXEC_DIR/bin/update_database --db-dir=$DB_DIR/ $META $COMPRESSION

Update 1: I was informed that it might be a memory issue, I'll double the ram configuration and try again.

Ram has been increased from 2GB to 6GB. and I'm testing it again.

Update 2: Still doesn't work with 6GB ram but it's improved, now I see a few new lines until it broke again. Increasing the ram to 12GB and testing again.

Update 3: It doesn't work even I increased the ram configuration to 12GB. 12GB has the same result as 6GB was. I gave up. Actually I think it might be the source file, did anyone has tried to import china map ever?

I figured. But who would have expected that importing a 400MB+ file can use more than 10GB of disk space? I increased the disk space to 200GB and Im gonna try again. I Hope this time it works.

Conclusion, if you are going to import china map, make sure you have 15GB+ free disk space and 4GB+ ram with your vm instance. I did it at last.

Also, some of the steps on the wiki telling you how to install overpass api was wrong.
http://wiki.openstreetmap.org/wiki/Overpass_API/Installation#Setting_up_the_Web_API

nohup ../src/bin/init_osm3s.sh $PLANET_FILE $DB_DIR $EXEC_DIR &
tail -f nohup.out

this is wrong, $PLANET_FILE was used as a folder and a file in different steps. here $PLANET_FILE is used as a file.

Better to reference another installation guide here to cross reference
http://overpass-api.de/full_installation.html


texts: Clustering features based on their proximity?

texts:
I have around 500 features (Forests) that I want to cluster based on their proximity, and then work with those clusters. I only have the ID, and shape area of those features.

How can I cluster them?


texts:
I have a dxf file and use this tool to convert it to shapefile.

Then I apply some affine transformations and get a shapefile that does not fit the shapefiles downloaded from osm (EPSG:4326).

I know that I have to reproject the shapefile using ogr2ogr (am I right?), but as I began without no clue of what was the original CRS, how can achieve the desired result of all districts in green of the picture on the right to fit as the magenta ones on the left one?

Thank you in advance !


texts: Arcpy SelectLayerByAttribute for loop with "ADD_TO_SELECTION"

texts:
I am a little bit confused with the "ADD_TO_SELECTION" possibility. I process a point layer. First, the idea was to process the points with the different GBAUPS attribute separately. And it worked like this:

periode_field = "GBAUPS"
valueList = []
rows = arcpy.SearchCursor(input_gws_layer)

for row in rows:
  valueList.append(row.getValue(periode_field))
uniqueSet = set(valueList)
uniqueList = list(uniqueSet)
uniqueList.sort()
print uniqueList

for period in uniqueList:

  period = int(period)
  bauperiod = arcpy.Select_analysis(input_gws_layer, "input_gws_layer_%i" % period, "GBAUPS = %i" % period)

  ...... and the processing ......

But I would like now to add continuously the points from the previous layers to the next one (first just one value of the attribute GBAUPS, then this one and the second one etc.)
I've tried like this, but the result is exactly the same as in the previous case.

for period in uniqueList:

  period = int(period)

  # selection
  arcpy.MakeFeatureLayer_management(input_gws_layer, "input_gws_layer_lyr")
  arcpy.SelectLayerByAttribute_management("input_gws_layer_lyr", "ADD_TO_SELECTION", "GBAUPS = %i" % period)
  arcpy.CopyFeatures_management("input_gws_layer_lyr", "input_gws_layer_lyr_%i" % period)
  ...... and the processing ......

I am not sure if I do something wrong in the SelectLayerByAttribute_management or in the CopyFeatures_management.


texts: How can I avoid creating a new cursor object repeatedly while counting features in a large SDE Database?

texts:
I would like to get a script running much faster than it is now. The script completes on my machine in about 33 minutes. It generates a report that counts features with specific attributes in one field.  What I think is generally the problem to be solved is that there are 7 features classes in 70+ datasets, with 4 possible entries in the field that I want counted separately.

Here is a very reduced and plain version of the code, it shows generally the structure of the nested cursor in some nested for loops.  I wonder if it is possible to NOT have the cursor creation be this nested??

state_plane_list = [...]
feature_class_list = ['PRCL_POL_',...]
clause = ["RESEARCH_STATUS = 'NULL'",...]

for state_plane in state_plane_list:
     for feature_class in feature_class_list:
          for rstatus in clause:
               cursor = arcpy.da.SearchCursor(feature_class + state_plane, fields, str(rstatus))
                       ***Counting code here - everything ends up in a dictionary and eventually gets written to a csv using DictWriter() ****

Basically, is there another...hopefully more efficient... way to do this?? I feel like creating 1600+ cursor objects may be way too dense a procedure.


texts: Adding custom constraint to a GeoDjango table?

texts:
I wish to add a custom constraint to a GeoDjango table.

This is the corresponding sql-statement:

ALTER TABLE my_table
   ADD CONSTRAINT enforce_valid_geom CHECK (ST_IsValid(geom));

Do I really have to inject the sql directly?

That seems like a frequently occurring constraint to me.


texts: Keep up to date a cpp-plugin(windows) for QGIS with less effort

texts:
We have two c++-plugins for QGIS and thought with the LTR release we are bit saver against the fast version-changes. Today I updated my QGis from 2.14.0 to 2.14.4 and had to recognise, that the plugins couldn't been loaded. Could be the reason therefore the new QGIS is built and deployed with a newer GDAL/OGR?

Or is there for every build of QGis requiered to have a seperate build for the plugin?

So I have to download and update the sources again and recreate the VS Project for every subversion of the LTR?

  EDIT:

Found out, need new gdal_i.lib and spatialite_i.lib to compile, so I changed it in the Linker-Properties of my VS-Project.
Then I need the libs qgis_core and qgis_gui. Downloaded them with the qgis-devel and linked it. But now I get the following Linker error:

  Fehler    1058    error LNK2001: Nicht aufgelöstes externes Symbol ""public: virtual bool __thiscall QgsSymbolLayerV2::writeDxf(class QgsDxfExport &,double,class QString const &,class QgsSymbolV2RenderContext *,class QgsFeature const *,class QPointF)const " (?writeDxf@QgsSymbolLayerV2@@UBE_NAAVQgsDxfExport@@NABVQString@@PAVQgsSymbolV2RenderContext@@PBVQgsFeature@@VQPointF@@@Z)". QGISFunc.obj    tb4qgisplugin

But in code I can't find any hint, never use sth like DXF. Any idea?


texts: How can I restart the loop or skip the error?

texts:
I set up a loop to extract cloudiness in multiple places, but if there is no image in one place in the loop, GEE will report an error.I want to loop to a place where there is no image and the loop can continue to the next place.Below is my code:
for(var i = 0; i < 10; i++)
{
  var id = ee.Feature(listall.get(i)).getInfo()['properties']['FieldID'];
  var nm = ee.Feature(listall.get(i)).getInfo()['properties']['FactoryNam'];
  var addre = ee.Feature(listall.get(i)).getInfo()['properties']['FactoryAdd'];
  var Lat = ee.Feature(listall.get(i)).getInfo()['properties']['lat'];
  var Lng = ee.Feature(listall.get(i)).getInfo()['properties']['Lng'];
  var x1 = ee.Feature(listall.get(i)).getInfo()['properties']['Lng1'];
  var y1 = ee.Feature(listall.get(i)).getInfo()['properties']['lat1'];
  var x2 = ee.Feature(listall.get(i)).getInfo()['properties']['Lng2'];
  var y2 = ee.Feature(listall.get(i)).getInfo()['properties']['lat2'];
  var x3 = ee.Feature(listall.get(i)).getInfo()['properties']['Lng3'];
  var y3 = ee.Feature(listall.get(i)).getInfo()['properties']['lat3'];
  var x4 = ee.Feature(listall.get(i)).getInfo()['properties']['Lng4'];
  var y4 = ee.Feature(listall.get(i)).getInfo()['properties']['lat4'];

  var firm = ee.Geometry.Polygon([ [x1, y1], [x2, y2], [x3, y3], [x4, y4] ]);
  var center = ee.Geometry.Point(x1, y1);
  var pol = center.buffer(5000, 500);
  var outer = pol.difference(firm);

  var sr = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");
  var t1 = ee.ImageCollection("LANDSAT/LC08/C01/T1");

  //卫星图像筛选
  var date_rg = ee.DateRange('2013-04-01', '2013-04-15');
  var SR = sr.filterDate(date_rg)
    .filterBounds(firm)
    .sort('CLOUD_COVER')
    .first(1)
    .clip(pol);

  var T1= t1.filterDate(date_rg)
      .filterBounds(firm)
      .sort('CLOUD_COVER')
      .first(1)
      .clip(pol);


texts: Code completion for matplotlib and numpy libraries in qgis script

texts:
I would like to know if it's possible to get code completion working for the matplotlib and numpy libraries in QGIS Script.

For example, the following code runs fine in script console of QGIS, but there is no code completion for neither libraries:


texts:
Possible Duplicate:
  Ordnance Survey grids - download or generate?

Does anyone know if I can download a grid of the tile locations for the 1:10k, 1:25k and 1:50k Ordnance Survey tile locations?

I know you can download the tfw/tab files for each tile, but I do not have the images. Is it possible to create this grid using tfw/tab files without the image files?


texts:
I was really curious about the new 3D-Map-View function of QGIS3. Now, finally, QGIS3 just released but I'm missing the 3D Map View! It should be at View -> New 3D Map View but it isn't.(https://anitagraser.com/2017/11/25/intro-to-qgis3-3d-view-with-viennese-building-data/)

Is it depending on my system configuration (Linux Mint 18.3 based on Ubuntu 16.04, Kernel 4.10)? Or have I overlooked a setting?


texts: Export shapefile with "non standard" characters to csv and read in excel

texts:
Macbook air mid 2011 running OSX 10.10.5
QGIS 2.8.2 - Wien
Excel 15.13.1 (O365)

When I export a shapefile containing non-ascii characters to a csv file using QGIS (setting encoding to utf-8) and importing to excel I end up with weird characters. What encoding (export and import) should I use?

for example:
Darién becomes DariÃ©n

A workaround is uploading the QGIS exported csv file with utf-8 encoding to google docs and download it as an xlsx file but there must be a more elegant solution

The way I export the shapefile is by using right mouse button on the layer and choose save as. In the dialog window I choose csv


texts: Incorrect result when using v.net.distance on grass76

texts:
As shown in the image below I need to connect the points in yellow to the respective nearest point in brown via the shortest path on the road network in blue.

I decided to use v.net.distance on GRASS GIS 7.6 for this purpose as the documentation for it describes exactly what I need to achieve. However, as visible in the image, the path that was generated (in red) takes a longer route instead of the more obvious shortest path as it is supposed to. Not sure what I'm doing wrong here.

The commands I used are as follows.

v.net roads points=bs op=connect thresh=10000 output=net1
v.net net1 points=end op=connect thresh=1000 output=net2 arc_layer=1 node_layer=3
v.net.distance input=net2 output=network flayer=2 to_layer=3


texts: What is geom_way column in osm2po generated data?

texts: Could anybody tell me the meaning of geom_way column generated by osm2po?


texts: If elif ArcPy If feature class equals then alter field / delete field

texts:
I have multiple feature classes in a geodatabase. I want to run a script that checks if a specific feature class exist within the geodatabase and if it does, modify specific fields found in that feature class.I want it to go through and find each feature class listed in the script and modify the fields. I have a beginning script that I exported from ModelBuilder to get a starting point, but I'm not sure how to put this into a loop.

# Import arcpy module
import arcpy

# Local variables:
FLANGE = "FLANGE"
FLANGE_Select = "\\\\\\My Documents\\ArcGIS\\Default.gdb\\FLANGE_Select"
FLANGE_Select__2_ = FLANGE_Select
FLANGE_Select__3_ = FLANGE_Select__2_
CASING = "CASING"
CASING_Select1 = "\\\\\\My Documents\\ArcGIS\\Default.gdb\\CASING_Select1"
CASING_Select1__2_ = CASING_Select1
CASING_Select1__3_ = CASING_Select1__2_
CASING_Select1__4_ = CASING_Select1__3_
CASING_Select1__5_ = CASING_Select1__4_
CASING_Select1__6_ = CASING_Select1__5_

# Process: Select
arcpy.Select_analysis(FLANGE, FLANGE_Select, "")

# Process: Alter Field (4)
arcpy.AlterField_management(FLANGE_Select, "WORK_ORDER", "Order", "New Order", "", "", "NON_NULLABLE", "false")

# Process: Alter Field (5)
arcpy.AlterField_management(FLANGE_Select__2_, "DATE_INSTALLED",    "Installed", "New Installed", "", "", "NON_NULLABLE", "false")

# Process: Select (2)
arcpy.Select_analysis(CASING, CASING_Select1, "")

# Process: Alter Field
arcpy.AlterField_management(CASING_Select1, "BEGIN_SERIES", "Series", "New Series", "", "", "NON_NULLABLE", "false")

# Process: Alter Field (2)
arcpy.AlterField_management(CASING_Select1__2_, "END_STATION", "Station", "New Station", "", "", "NON_NULLABLE", "false")

# Process: Alter Field (3)
arcpy.AlterField_management(CASING_Select1__3_, "EVENT_GUID", "Guid", "New Guid", "", "", "NON_NULLABLE", "false")

# Process: Delete Field
arcpy.DeleteField_management(CASING_Select1__4_, "USER")

# Process: Delete Field (2)
arcpy.DeleteField_management(CASING_Select1__5_, "ID")'


texts:
I'd like to calculate the median width of a polygon in PostGIS.

The problem is almost identical to the one described here, which has an excellent answer. My question is how to calculate the median, not the mean.

The polygons represent roads, and due to irregularities, the median is a more robust measure than the mean. I have the road centerlines if that would help, although due to positional inaccuracies they are not always precisely in the center.

--UPDATE--

Here's an example of where the mean would be different from the median - a cul-de-sac where the mean would be wider than the median due to the turning circle at the end. The yellow polygons are the street rights of way.


texts: OpenLayers.Bounds transform not working as expected

texts:
I have problem with transforming bounds. It is not working for the first time, later it works just fine. At first i thought that something was still loading butsearching in firebug shows that nothing gets loaded and i think everything is provided for it to work:

Here is the code:

        Proj4js.defs["EPSG:2180"] = "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs";

        var proj2180 = new OpenLayers.Projection("EPSG:2180");
        var proj4326 = new OpenLayers.Projection("EPSG:4326");

        var bounds = OpenLayers.Bounds.fromArray(record.get("llbbox"));
        bounds.transform(proj4326, proj2180);
        console.log(bounds);

This code is placed in function. When i run it for the first time i get following console output:
13.897781645001,48.960416159343,24.621296204369,54.83756633314

when second time i get:
143457.09623933677,137520.1675486639,876176.7241084952,786832.4073041743


texts: Calculate area of two intersecting rasters

texts:
I have two large rasters: One defining zones (private properties, with unique value for each property) and one of land cover change (deforestation data, with values for the year it has been deforested).

What I want is to calculate the amount of deforestation for each year on each property. I will later categorize the properties by property size, but I do not want to do it at this stage. And I would like to be able to plot the results spatially as well.

Also, I am interest on the size of the deforested patches, i.e., adjoining raster cells deforested on the same year. Is there a way to calculate this area without transforming into polygons?

I know that using polygons I could use intersect and it would give me exactly what I'm looking for. But these are VERY large rasters, it is simply not possible to use polygons for this. I looked into Zonal Statistics but as values for deforestation are categorical, that doesn't solve it. I'm currently trying to fix a bug while doing Tabulate Area (backgoung processing crashes), but still I wouldn't be able to plot the results.

Is there any different approach for this? I feel it's a simple query that I just haven't found the right tool(s) yet...
I'm using ArGis 10.2.2


texts: How to transform a DEM to a cartesian coordinate system?

texts:
I am using the script r.inund.fluv. I run the script with all the required parameters and it produced a lot or errors.

I emailed the author of the script and he said that

"I think you are not using a cartesian coordinate system... ("The DTM resolution is 0:00:03 x 0:00:03 meters") and the code work with a cartesian grid. "

I'm using a raster map from SRTM for the elevation. Its DTM resolution is 0:00:03 x 0:00:03 meters. As what is required for the script to work, the author said that I should convert it to a cartesian coordinate system. How is it done for a raster SRTM map?


texts: Refreshing web appbuilder widget during development

texts:
I'm beginning to learn widget creating in Esri's Web Appbuilder for ArcGIS, but I'm having some problems getting the code to display when I make changes.

Coming from a Visual Studio-based web/Windows development background, I'm used to being able to make a few minor changes and then hit the debug button to see the changes in action.  With the widget development, I'll make changes to my code (widget.js or widget.html) but then I'm unable to view the changes in the Web App (it always shows the last visible version).

Simply refreshing the page does not work, nor does clearing the cache or other web data.  I've even tried restarting the node server and Web Appbuilder itself, but no go.  The only thing that has worked so far is to delete the app I put the widget into, then create a new app and load the widget into that one.  This seems to work most times, but it is tedious and time consuming.  I got the idea from https://geonet.esri.com/thread/120241, but that is from the beta of Web Appbuilder nearly 2 years ago.  Is there a better method for regularly refreshing a widget during development?


texts: Using Python Parser of Label Expression to only show text after particular symbol?

texts:
I have a bunch of point data that needs labeling. However the attribute used for labeling is more complex than is needed in the map.

Example of the attribute field for labeling:

15/8-9a,
001/15-2
8/19-13

What we need to see labeled is any number and/or letter found after the "-" symbol.
So from the above example: 9a, 2, 13 will be displayed on the map.

I know this can be done with slicing but since the number of characters after "-" is variable I'm struggling a little with the expression.

I do not want to add a new field so am exploring writing a label expression.

Does anyone know what the expression might look like to handle this?


texts:
I have 5 different polygon vectors:

All have over 20.000 features;
All have multiple attributes;
All overlap each other at least more than 80%, but are shaped differently;
All display information of large areas, about 1/3 of the North Sea;
All have an attribute called "benefit" which ranges between 0-1 (I normalized the data);

My mission: Using the attribute "benefit" I want to find optimal locations for a specific activity. The sum of attributes for 5 different vector files, would indicate where I can find such optimal sites. The higher the better.

For example, let's say I choose two different geographic locations close to each other somewhere off the eastern coast of England. For each location I would retrieve the relevant value of the attribute "benefit".

For the attribute "benefit" I would get:

Location 1
layer 1: 0.5
layer 2: 0.3
layer 3: 0.6
layer 4: 0.8
layer 5: 0.9

Location 2
layer 1: 0.1
layer 2: 0.3
layer 3: 0.2
layer 4: 0.1
layer 5: 0.4

Location 1 would be better compared to location 2.
Sum layer 1: 3.1 vs Sum layer 2: 1.1

THE QUESTION: What would be the best methdos to do this for all layers across the all extends?

At my disposal I have: QGIS, Python and pgAdmin

1) QGIS would be the least technical but I'm not even sure it's possible. With such large vector files I think it would crash. I could use the python console maybe...?

2) pgAdmin might be quicker, but for each feature I would need to check if it overlaps any of the other features...and with 5 vector files each containing more than 20,000 features...(some even 100,000)...I'm not sure if this is possible.

3) Python might also allow for a quick analyses of the data. But I'm not sure which packages to use.


texts: QGIS - Is it possible to display a raster label (not band value) using identify tool or any other tool?

texts:
This question was asked 4 years ago, but I hope that new QGIS features may have produced an anwswer.

In the new QGIS version, are there new raster features  similar to the vector features "action button", label menu or Qt Designer form, that may allow  to display a raster label (not band value) ?


texts: Appending 10,000 tables using ArcPy fails to complete

texts:
I have been attempting to append 10,000 tables using the following arcpy code (from an answer to Appending 10,000 tables using ArcPy?):

WD = "C:/"
arcpy.env.workspace = WD

for table in arcpy.ListTables():
    if not table.startswith('test_4'):
        tableData = list(arcpy.da.SearchCursor(table, "*"))
        cursor = arcpy.da.InsertCursor("test_4", "*")
        for data in tableData:
        cursor.insertRow(data)

This code works perfectly for 5 or 500 tables. When I try to append 10,000 tables, the program runs for 5 minutes and then stops running. It does not give an error message, however. The program seems to still be running, although it is in fact not.

What could be causing this?


texts: Road Type Symbology Query (non technical)

texts:
Does anyone know of a good resource/resources for road type symbologies?

I am looking to update the symbology of a Road class layer which has 8 categories. The Road schedule is for Ireland. So a standardised symbology for the EU/EU country would be ideal. I know road types would differ in differ countries but some ideas would be great!

I am currently using symbology from Ordnance Survey Ireland by as you can see the Local Primary-Tertiary symbologies are pretty much indistinguishable.


texts:
I want to integrate QGIS 2.12.3 or 2.13 with GRASS 7, but I still haven't been able to do so.

Can someone please tell me how can I do this with a stable QGIS version.


texts: I have maps generated using ArcGIS (.tfw, .fig, .prj, .sbn, .shp and .dbf files).  I need to convert them to ADRG format.  Can ArcGIS do that?  If so, how?


texts:
I am trying to mosaic all the data from VIIRS of a day into a single gtiff.

I've already done it with MODIS data using GDAL. The problem it seems is that the coordinates of each h5 files doesn't seem to be in the same file as the data. Which gives this as an error when I try a gdalwarp:

  ERROR 1: Unable to compute a transformation between pixel/line and
  georeferenced coordinates for

  HDF5:SVI01_npp_d20140812_t0007456_e0013260_b14451_c20140812061326517212_noaa_ops.h5://All_Data/VIIRS-I1-SDR_All/Radiance.

  There is no affine transformation and no GCPs.

Do you have an idea what could be done to get the localisation data and to merge them together ?


texts: How to emulate this grayscale basemap look+feel in Google Maps

texts:
I love the design on the Development Seed homepage, where they've used a grayscale basemap with bold 3D-style graphic markers.

It appears that they've used the ModestMaps API to create a tiled basemap at various resolutions, with the text and graphics burned into the basemap.

I'd like to emulate this look+feel (for a personal, just-for-fun project) using freely available tools such as the Google Maps API. I need the ability to change the marker overlays rapidly, so they would need to draw dynamically over the basemap.

Is it possible to achieve a similar grayscale, terrain-based backdrop using Google Maps?

The Google Maps API Styling documentation shows how to reduce the saturation of the topo map, but I can't see an option to do it with the Terrain map.

I think I can emulate the bold graphics using the complex Icons sample as a starting point but would welcome any suggestions.

Thanks


texts: I am creating a Heat Map and need Conservation Internationals's Biodiversity Hotspot Data in GIS form. I have tried to find the data online, but have since been unable. Where would I be able to find this data?


texts: Choosing towns with largest area with exact amount of inhabitants

texts:
I have polygon layer of Villages and I have to choose an area of Villages as large as possible that have (together) the same amount of inhabitants as the largest one.

Basically, I need to create a layer of Villages with area large as possible and with same amount of inhabitants as the Capital.

Now I see that i need to use "knapsack problem" but I'm not skilled in Python. Is there a method in ArcMap (model or something) that can solve this?


texts: Reference system with all grid lines being great circles

texts:
I was thinking of a problem where it might be easier if instead of latitude (and parallels) we would have a system with an additional east/west pole and great circles through them would provide a second coordinate to go along with longitude.
Was thinking surely this has been done before, but I'm missing a name for such a reference system.
Is this known/used and by what name so I can read more?
edit: a sketch as requested:


texts: How to turn an IField from Not Editable to Editable?

texts:
so i have turned my ITable into an editable DataGridview from this post
link.

Now i want to set some of the columns to be non editable. So what I can do is initially set my IFieldEdit.Editable_2 to false, and the column would be read only.

However, when I tried to do that and actually populate values into my field initially, it throws me exception for not editable.

So I was wondering if anyone know if i can temporary set the field to editable and then after setValue , set the editable to false again.

I have this but doesn't work:

bool fieldEditable = featureClass.Fields.get_Field(x).Editable;

IFieldEdit fieldEdit = (IFieldEdit)featureClass.Fields.get_Field(x);
fieldEdit.Editable_2 = true;

featureBuffer.set_Value(x, cegObject.GetPropertyValue(featureClass.Fields.get_Field(x).Name));

fieldEdit.Editable_2 = fieldEditable

However, the setting to true, does not work. Is there a way around this? Or is editable totally useless cause once i make a field editable, i can't even populate a value in? But i'm sure i'm just not doing it right, please help.

Thanks and Regards,
Kevin


texts: Converting an annual Image Collection into a multiband image in Google Earth Engine

texts:
I have an image collection with 35 years of imagery. Each image in the collection represents a measure of annual NDVI for each year from 1985 to 2019. I want to take this image collection and create a single image where there are 35 bands, each band representing NDVI_1985... to NDVI_2019.
I think it would something like this:
var ndvi_collection = ndvi_coll.select(['nd_p95_rescaled','time'],['NDVI','time']);
print(ndvi_collection, "ndvi collection");

var multiband = ndvi_collection.iterate(function(image, result) {
  return ee.Image(result).addBands(image.select(['NDVI']));
    }, empty);

The issue is that I want to carry over information about year in the band names or retain the 'time' band for each year as well.
Here is the link to the entire GEE script: https://code.earthengine.google.com/a7fcc54f529ca1a43743f90b603a85b7


texts: How to display/render OGC WFS web service on Google Maps using Google Map JavaScript API v3

texts:
I want to display publicly available OGC web services on Google Maps. I figured out the solution for WMS but no clear solution for WFS.

Sample WFS: http://services.ga.gov.au/geoserver/boreholes/ows?service=WFS&request=GetCapabilities&


texts: Extract percentile info from multi-polygon image

texts:
I'm trying to find the lowest 20th percentile figure and corresponding area across different geographies defined by an image.
The image is obtained from the multiplication of 2 images, which has produced many disparate polygons in brown (soil taxonomy and slope). See below.
var img1a = soiltaxa.eq(389).and(slope0206Mask);
Map.addLayer(img1a.updateMask(img1a),{palette:'brown'},'1a');

Separately, I created a single band for soil carbon below in black.
var SOC = socc.multiply(bd)
  .multiply(5*0.1*10000/(250*250))
  .clip(region);
var imgSOC30 = ee.Image(SOC.select('b0','b10','b30'));
var SOC30 = imgSOC30.reduce('sum');
Map.addLayer(SOC30,{min:12,max:242},'SOC30',0);

I want to find the lowest 20th percentile and corresponding area of the single band represented in the bottom picture in monochrome (soil carbon) for just the brown geographies in the top picture.
The issue seems to be that the brown polygons in the top picture are an Image. How do I use Google Earth Engine to calculate the the lowest 20th percentile soil carbon and corresponding area for just the brown polygons?
The code below works if I substitute "img1a" for a geometry, but when I use "img1a", I get an error message saying:
Dictionary (Error)
Image.reduceRegion, argument 'geometry': Invalid type.
Expected type: Geometry.
Actual type: Image<[grtgroup]>.
var SOC30percgeo = SOC30.reduceRegion({
  reducer: ee.Reducer.percentile([20,50,80]),
  geometry: img1a,
  scale: 250, // Resolution in meters.
  maxPixels: 1e9});
print('SOC30 percentiles for geo',SOC30percgeo);

// Get lowest 20th percentile SOC30 for geo.
var DegradedSOC30geo = SOC30.clip(img1a).lte(ee.Number(SOC30percgeo.get('sum_p20')));
Map.addLayer(DegradedSOC30geo.updateMask(DegradedSOC30geo),{palette:'red'},'Degraded SOC30 geo');

// Get pixel area for geo degraded and convert to ha.
var DegradedSOC30geoArea = DegradedSOC30geo.multiply(ee.Image.pixelArea()).divide(10000);

// Sum pixels for geo degraded.
var DegradedSOC30geoSize = DegradedSOC30geoArea.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: img1a,
  scale: 250, // Resolution in meters
  maxPixels: 1e9
});
print('degraded area in geo: ',DegradedSOC30geoSize,'ha');


texts:
I have a LineString table in PostGIS, which represents wiring in streets. Where the geometry column represents the wiring lines; firstpoint is the first geometry point (cable segment); lastpoint the last geometry point; id_startpoint and id_enpoint are identifiers of the firstpoint and lastpoint columns.

In OpenLayers v5 it is shown as follows:

The point is that I want to detect the rings (cycles) in the structure of the network. As shown below:

My solution was to use the DFS algorithm in javascript to detect cycles in undirected graphs, which works well until I give the 5327 registers (too much data), however with 200 registers it works fine. For this reason, I want to use some PostGIS function to find the cycles. You could recommend some function to me since there is ST_NumInteriorRings and ST_ExteriorRing but they are for polygons, I need a function for LineString.


texts: QGIS: Clip a shapefile based on values in a raster

texts:
I would like to clip a shapefile based on the values of a raster file beneath it, using QGIS.

The specific context is a shapefile of Local Government Areas overlaid onto a population density raster. I would like to clip the shapefile based on the values of the raster file, e.g. clip the shape Local Government Area to only include areas with a population density above 50 people per square kilometre.

I have found resources on how to do the reverse (clip a raster from a shapefile) but not to clip a shapefile based on raster values.


texts: Positioning of graticule labels in QGIS 2.14.1

texts:
I have a map grid (line, shp) that was created with MMQGIS.
The X and Y coordinates of the lines are stored in the the attribute table, which makes labelling the lines an easy process.

My problem is with the location of the labels.
Currently the labels are in the default position, central. This poses an issue when I want to print a gridded map as the labels are very prominent and obscure information in the center of the map.

I wish to align the labels of horizontal lines to the left margin of the page and labels of vertical lines to the top of the page.

I have seen an identical question asked previously here SHP file Left / Right Label Alignment issue in QGIS 2.0.1-Dufour
I have attempted to use the same solution with no success.
I am using version 2.14.1 and the coding required to position the labels at the start or end of the line must be different, or I am using the code incorrectly.

I have attached an image as an example of the default. I wish to move the labels to the left for the horizontal lines, and to the top for the vertical lines.
If it is possible to offset the labels from the edge of the page that would also be great.

If any more information is needed please let me know.

Image 01: Centered labels

Edit 01, for clarification:

I have managed to align the labels to the edge of the lines, although this is at a fixed position (see image).

The problem is that if I zoom in on a specific area I cannot see the labels.
I would like the labels to follow the screen boundaries if this is at all possible. (It may not be currently)
By doing this, it should hopefully automatically add the labels to the edge of a map created in the print composer.

Image 02: Zoomed out image, labels visible.

Image 03: Zoomed in image, labels in fixed location, therefore not visible.


texts: How to change overlays checkbox feature in layerlist to radio button

texts:
How to change the default check box button of Overlays in Layerlist to radio buttons. My client doesn't want his users to see more than one layer on the map.

Is there any way that we can prevent the possibility to check multiple layers in overlay layers in the LayerList?


texts:
I'm trying to view some data held in a SQL Server 2008 database using MapServer. I can run the GetCapabilities requests without any errors, but when i try and run a GetFeature I get the following error:

<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.2.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
<ServiceException code="NoApplicableCode" locator="mapserv">
msWFSGetFeature(): WFS server error. ms_error->code not found msMSSQL2008LayerOpen():
Query error. Couldnt make connection to MS SQL Server 2008 with connect string 'Server=xxx;user id=yyy;password=zzz;database=db;trusted_connection=yes'.
Error reported was '[Microsoft][ODBC SQL Server Driver][SQL Server]Cannot open database "db" requested by the login. The login failed.'.
This error occured when trying to make a connection to the specified SQL server.  Most commonly this is caused by
(1) incorrect connection string
(2) you didnt specify a 'user id=...' in your connection string
(3) SQL server isnt running
(4) TCPIP not enabled for SQL Client or server <br>
</ServiceException>
</ServiceExceptionReport>

When I use the same connection string in a VB.Net aspx page, I can get into the database and table I'm trying to access with MapServer, so I think the connection string is ok. I've enabled the TCP/IP protocols in the SQL Server Configuration Manager. I think this is a MapServer issue in that it can't get to my SQL 2008 database in the same way as my aspx page can.

I've also tried setting the SQL 2008 layer connection details using the "plugin" connection type and also as an OGR connection in my mapfile but without any luck.

Can anyone give me any pointers as to what to try next please?


texts:
i have configured a tile server by following the switch2osm

in the last step for starting the render automatically on startup
if i start the render as below vars.sh is missing

$ sudo /etc/init.d/renderd start
root's password:
/etc/init.d/renderd: line 33: /lib/init/vars.sh: No such file or directory

What is have missed in this?


texts: ArcMap 10 - Viewshed Analysis - Input Type

texts: Hopefully someone can tell me if I am missing something here. I am hoping to built a TIN in ArcScene and run some view shed analysis on it, simple enough but I need to add building heights to the TIN, this was done through an ArcGIS Script. The problem is I can't seem to get the TIN as the input surface using ArcMap 10 viewshed anlysis in the 3D or Spatial Analyst. Has ESRI removed this parameter completely from version 10?


texts: Visualization of dynamic data on top of OpenStreetMap tiles?

texts:
I want to show dynamic data (say .2 sec update) from a postGIS database displayed on top of OpenStreetMap data.

I am very new to GIS so not sure how to go about this.

I am thinking something like:

Prerender OSM tiles + serve with GeoServer (or MapServer).
Serve the dynamic data separate from GEOServer. This could just be lat long coords.
Use OpenLayers in broswer to place the dynamic data on top of the tiles.

Is it possible to do what I am thinking here and if so, does this seem reasonable?
Would it be better to render everything server side and send back result? Or should I not prerender tiles?


texts: MapProxy dosen't display certain uncached tiles

texts:
I face an unnormal behavior of MapProxy on client.

The way of spatial data looks like: |PostgreSQL| --> |GeoServer|wms request| --> |seeding tiles|MapProxy|wmts request| -->|Fronend(OL6)|

My seed.yaml looks like:

seeds:
  myseed1:
    caches: [ps_osm_cache]
    levels: [0,1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
    refresh_before:
        time: 2019-12-04T12:35:00

But some tiles on 13-14 levels:

12 (all is ok):

13 is not (water styles are dead):

14 bad again (but some new correct tiles):

and from 15 to 19 all is OK!

What is wrong? It isn't mapproxy cache of 13-19, cause I tryied remove all tiles, and did seeding again... result hadn't achived...
And It isn't 'no data' problem, cause wms from geoserver and PostgreSQL info show that all is ok.. I really don't understand what the reason..


texts: QGIS: Saving TIF as Geopackage raster loses NoData cells

texts:
Background:

At 3.8.3, I right-click a TIF DEM (single-band, 32-bit floating point) and Export > Save As to a Geopackage raster.

If I use all default settings in the Export > Save As window, the target raster is identical to the source TIF.  I am happy.

Problem:

However, if I select the Pyramids option in the Export > Save As window (with External GTiff .ovr, cubic resampling), the NoData values in the input TIF  are incorrectly set to the value -999999.1 in the output Geopackage.  I am sad.

I also tried to force NoData output by checking the NoData option (screenshot below) and inserting a range that would encompass the -999999.1 range.  No difference.

Workaround:

If I Export without the Pyramids option, then run the Raster > Miscellaneous > Build Overviews tool, the resulting NoData cells correctly remain as NoData.  I am happy.

Question:

Is the NoData -> -999999.1 output normal, or (as I suspect) is this a bug in the Export > Save As window?


texts:
I have a browser-based work order management system (Maximo 7.6.1.1).

I want a user to be able to execute a URL in Maximo that would open a web map and filter a layer.

For example, tell the web map to filter a layer based on a list of IDs:

wonum=WO1778,WO1772,WO1763

How can I trigger the Portal web map to filter a layer via a URL?


texts:
I have some metadata saved using the ArcCatalog metadata editor.  The abstract element allows you to use html formatting.  When you export the metadata to xml the formatting is stripped out and the text is returned as a single line of text with newlines removed.

Is there a way to programmtically access these elements from Python and retain formatting and newlines.

Ideally this would be as an XML document.


texts: Overlaying city and borough names with political boundaries in QGIS?

texts:
This question is related to Fine-Grained Political Modern Map of Italy.

The following images are from  @Weinfried's answer. I was able to get the ITA_adm0.kmz, ITA_adm1.kmz, ITA_adm2.kmz, but not ITA_adm3.kmz, which for some reason does not open (I extracted these prior to importing them into QGIS). However, I don't know how to view the city names from Google Maps or other source in QGIS. Can someone please help me out?


texts: QGIS 2.10.1 Clipping Algorithm Results In Empty Layer

texts:
I apologize for the redundancy of this issue, but I have explored all the information I could find over the past couple days and have yet to find a solution to what seems to be a very simple task.

I am attempting to clip a larger landscape layer to a smaller clip layer in QGIS.  The resulting layer contains the attribute column of the source layer, but no features.  I have run the geometry validity test on both layers with zero errors returned.  I have ensured both layers are of the identical CRS.  I have experimented with expanding the extents of the clip layer to overlap those of the landscape layer.  Everything I have tried results in the same empty layer.  The same issue exists when attempting to perform the intersection operation.  After experimenting with other layers, I have determined the issue to be with the landscape layer and none other.  HERE is a link to the project files so anyone may explore them further.

EDIT: I took another look at the OGR error log and it repeats the error: "Feature geometry not imported (OGR error: Attempt to write non-polygon (GEOMETRYCOLLECTION) geometry to POLYGON type shapefile.)" for what I presume is every feature to be created in the clip output, but is not due to whatever this error means.  Anyone care to shed some light on this?


texts: What error am I doing when rendering OpenStreetMap roads and buildings in Minecraft?

texts:
I'm trying to create a Minecraft map of the Ernst Reuter Platz in Berlin. For this purpose, I'm using OpenStreetMap data.

This is what the place looks like on a 2D map:

When I render the roads and the buildings into a Minecraft map, I'm doing something wrong. Here you can find a zoomable preview of the map.

Gray blocks represent streets (from highways to footways) and reddish-brown ones - contour of the buildings.

As you can see, there is some error in the renderer. For example at a certain location there is a building, where the lines seem to be drawn in the wrong order.

Is there an algorithm, which I can use in order to sort the points of the OSM polygons so that they result in a nice shape (like a real building) ? If yes, what is it called?

Update 1 (30.09.2015 08:38 MSK): One possible approach is to use convex hull algorithm.

Update 2: Sample definition of a building

<way id="103989044" visible="true" version="4" changeset="34001855" timestamp="2015-09-13T15:17:07Z" user="species" uid="321042">
    <nd ref="1200207407"/>
    <nd ref="1200207412"/>
    <nd ref="1772686869"/>
    <nd ref="1119746578"/>
    <nd ref="1119746577"/>
    <nd ref="1772676000"/>
    <nd ref="1200207407"/>
    <tag k="addr:city" v="Berlin"/>
    <tag k="addr:country" v="DE"/>
    <tag k="building" v="yes"/>
    <tag k="layer" v="1"/>
    <tag k="min_level" v="1"/>
    <tag k="source" v="Bing"/>
</way>


texts: PostGIS & Python SQL Query formatting: Create Geometry only for entries that are not NULL

texts:
I think this is rather a formatting issue. I have a database with a geometry Column that sometimes contains points. Normally, it will be None/Null.

When I format the SQL query in psychopg2 and python (for inserting entries), it looks like this:

insert_sql = '''
        INSERT INTO "post" (post_guid, post_latlng)
        VALUES (%s,ST_GeomFromText(%s,4326))
        ON CONFLICT (post_guid)
        DO UPDATE SET
            post_latlng = COALESCE(EXCLUDED.post_latlng, "post".post_latlng),
        '''
cursor.execute(insert_sql,(post_guid,photo_latlng))

However, this does not work if photo_latlng = None
I tried:

if geoinformation_available:
    photo_latlng = "POINT(%s %s)" % (photo_longitude,photo_latitude) #WKTconversion
else:
    photo_latlng = None
if not photo_latlng == None:
    photo_latlng = "ST_GeomFromText(" + photo_latlng + ",4326)"
insert_sql = '''
        INSERT INTO "post" (post_guid, post_latlng)
        VALUES (%s,%s)
        ON CONFLICT (post_guid)
        DO UPDATE SET
            post_latlng = COALESCE(EXCLUDED.post_latlng,"post".post_latlng)
        '''
cursor.execute(insert_sql,(post_guid,photo_latlng))

.. but it gives me a psycopg2.InternalError: parse error - invalid geometry
(HINT:  "ST" <-- parse error at position 2 within geometry). I also doubt that this will work with COALESCE -> only fill row if entry is not Null.


texts: ST_MapAlgebra and combination of four rasters

texts:
I am trying to create combination of four rasters, but I am lost with ST_MapAlgebra callback function. Official documentation (http://postgis.net/docs/RT_ST_MapAlgebra.html) gives an example, that return only 0 but does not read any values from input.

How should I create, for simplicity, callback function that will sum together four input rasters?


texts:
I have spatial function in my Postgres/PostGIS DB that when sent a set of coordinates (Long Lat) returns single value. Whilst I'm not a web developer, i'd like to serve the function up as a web service?

GeoServer and Apache Tomcat are currently installed on the same server and I wanted to understand if there is way to use the existing infrastructure to call the function?

Any advice/help greatly appreciated.


texts: Splitting large shapefile into numerous smaller files based on number of objects in each smaller file using?

texts:
I have a shapefile containing approximately 600 000 point-objects.

Unfortunately, for a couple of reasons I need to split this file into a larger number of shapefiles each containing approximately 200 objects.

Is anyone aware of any tools/scripts that allow me to automate this process?

I have tried to automate the process using e.g. the Split-tool in ArcGIS, however since I am more of an archaeologist than a programmer (and know next to no Python), I haven't been able to produce a script that does this.

I am guessing that the best way to do this would be to use a counter which increments with +1 for each object passed into a new shapefile, and when that counter reaches 200 it resets and a new shapefile is created into which the subsequent 200 objects are copied.

However, I haven't been able to write a script/create a workspace that does this.

So, does anyone have any suggestions of scripts/programs that are able to do this?

I am using FME Workbench.


texts: All nodes are duplicate in polygonize result

texts:
I have a strange error with a large complex polygon shapefile (representing a hypothetical historic flood level of a lake) that I cannot resolve.  Using QGIS 2.18.14, the polygon was derived as follows...

Reclassify a digital elevation model (DEM) into two classes- one class above the desired flood level, one below the desired flood level.
Polygonize the class below the desired flood level (Raster--> Conversion--> Polygonize).
Subsequently, I edited the geometry of the resulting polygon to merge small islands with the larger lake level polygon.  I also adjusted some boundaries manually using the edit nodes tool.

Now, when I open the shapefile, toggle editing, and use the nodes tool, validation starts as usual, except that the entire polygon consists of duplicate nodes.  Because the polygon is large and complex (~35,000 nodes), it takes some time to validate and is really inefficient to manually inspect.  I've tried the following remedies:

Use the GRASS v.clean tool
Use the MMQGIS Modify--> Remove Duplicate Geometries Tool
Import the shapefile into SAGA and convert the polygon to edges + nodes, then subsequently remove duplicate points with the Shapes - Points--> Remove Duplicate Points tool
Export only the edges from the previous step, saving them as a polygon
shapefile

None of these remedies seem to work.  I've not been able to find this same problem posted elsewhere online. I would ordinarily give up at this point and try to re-create the flood level polygon to better document what may have caused it to go wrong along the way (I suspect something happened during my manual editing, but I can't imagine what).  I put several hours into step 3 above (manual editing), and so it would be great if there was a way to solve this issue.

I included a link to download the shapefile below if anyone wishes to play with it.

p.s. - you may be tempted to criticize why I am even dealing with this type of analysis by converting an artificially flooded DEM into a polygon... it is a client's request to have the result as a shapefile.

Lake polygon download


texts: pgrouting crashes and restarts the postgresql server when searching for endpoints in 2 disjoints graphs

texts:
I have

>apt-show-versions postgresql-9.1-pgrouting
postgresql-9.1-pgrouting/precise uptodate 2.0.0-release-ppa1~precise2

Whenever I do a search http://docs.pgrouting.org/dev/src/astar/doc/index.html between two end points in disjoint graphs, I get  the error error sending I/O to backend and my postgresql server crashes and restarts. Why doesnt this simply return a no route found?


texts: Removal of Qgis 2.2 after 2.4 install, while keeping WFS and WMS features

texts:
Can i safely uninstall Qgis 2.2 after I've installed the 2.4 version and keep my WFS and WMS features?

I'd rather be certain, so I don't have to setup the whole thing again.


texts: Displaying nearby locations as 1 point (QGIS 2.8 Wien)

texts:
I am using an excel file to store all my locations and their attributes.
This excel file contains coordinates of the same location measured on different dates. This means that the coordinates are not exactly the same for every locality visit.
Is there a way that I can display these nearby locations as 1 single point?

To explain with an example: If I visit the zoo every weekend and save a coordinate of the zoo every time, the saved coordinates will not be exactly the same every time... Each weekend I add in my excell file what animals I have seen. Is it possible to show the zoo as one point (or one central coordinate) on my map and still show in the attribute list what animals I have seen during ALL of my trips and which animals I have seen at what date ?


texts: Looking for the census blocks that make up the most recent U.S. Congressional Districts

texts:
I am hoping someone can help me solve this or point me in a direction. I need to know what census blocks make up the most recent Congressional Districts (115th). I have found a file that the US Census Bureau has which is the census blocks as they relate to older congressional districts (2013 I think), but not the most recent ones. The closest I can get is census tracts inside of congressional districts and I can work down from there. My understanding is that tracts can be cut by a congressional district (Which census geography boundaries do congressional districts preserve?), and I am not sure how to solve that.

I have tried an intersect, but the boundaries do not seem to work out quite correctly and I get some weird things included/excluded. A geoid relationship based text file would be ideal.


texts: Installing postgis on postgresql 9.1 not working

texts:
I'm trying to install Postgis 2.0 on a perfectly functionig postgresql 9.1.10 database.

I installed the package via:

richard@srv1:/$ sudo apt-get install postgresql-9.1-postgis postgis
Reading package lists... Done
Building dependency tree
Reading state information... Done
postgis is already the newest version.
postgresql-9.1-postgis is already the newest version.
0 upgraded, 0 newly installed, 0 to remove and 2 not upgraded.

With out any errors.
Then I try to add postgis as an extension to an existing (empty) database, with:

 dev=# CREATE EXTENSION postgis;
 ERROR:  could not open extension control file "/usr/share/postgresql/9.1/extension/postgis.control": No such file or directory

If I look at the location the file does not exist.

What am I missing? I can't figure it out. If I google I only get ancient results....

Where can I find the version of the installed postgis? As I'm not able to do a SELECT PostGIS_version();

Environment: Postgresql 9.1.10, Ubuntu 12.04

Update

I'm willing to install postgis2.0; apparently I have 1.5 installed. After purging that version I added the repository

deb http://apt.postgresql.org/pub/repos/apt/ precise-pgdg main

After trying to install the package via

sudo apt-get update
sudo apt-get install postgis

it ends with:

Errors were encountered while processing:
/var/cache/apt/archives/libgdal1_1.9.0-3.1~pgdg12.4+1_amd64.deb
E: Sub-process /usr/bin/dpkg returned an error code (1)

Update 2

It appears because of the naming of the packages, it installed the wrong architecture type package.
But I ended up using this website and these steps after making sure I had gdal>= 1.9.0 and GEOS >=3.3.3:

wget http://postgis.refractions.net/download/postgis-2.0.0.tar.gz
tar xfvz postgis-2.0.0.tar.gz
cd postgis-2.0.0
./configure --with-gui

make
sudo make install
sudo ldconfig
sudo make comments-install

sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/shp2pgsql
sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/pgsql2shp
sudo ln -sf /usr/share/postgresql-common/pg_wrapper /usr/local/bin/raster2pgsql

TADA

durham=#  select PostGIS_full_version();
NOTICE:  Function postgis_topology_scripts_installed() not found. Is topology support enabled and topology.sql installed?
                                                                     postgis_full_version
--------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER
(1 row)


texts: Where is the official website of the JTS project?

texts:
I'm using the JTS API, published at http://www.vividsolutions.com/jts/JTSHome.htm
However, there is another JTS API, published at http://sourceforge.net/projects/jts-topo-suite/

The former is version 1.8, from 2006.
The latter is version 1.13, from 2013.

I don't understand. What is the right version of JTS to use?
If the right version is at www.vividsolutions.com, then is this project alive? There are no new features there since 2006.


texts:
I'm having trouble optimizing this query, sometimes it will take around 10s to execute. I have indexed the columns samples.geom and jbb2016.geom_3021 so that shouldn't be a problem.
I know this is surely a common question, but I have not found a solution to this slow query in my searchings.

Here's some info about my PostGIS installation: POSTGIS="2.3.0 r15146"

PostgreSQL version: PostgreSQL 9.6.1 on x86_64-apple-darwin14.5.0

The samples table contains about 400k rows the geometry is points, the jbb2016 contains about 1.25m rows the geometry there is some kind of multipolygon.

SELECT
year, customer_id, COUNT(customer_id) AS samples,
(SELECT name || ' ' || last_name FROM customers WHERE id = customer_id) as name
FROM samples
WHERE
ST_Intersects(
    (SELECT geom_3021 FROM jbb2016 WHERE blockid = 62023246713), samples.geom
)
GROUP BY year, customer_id
ORDER BY year DESC;


texts: Adding points in Empty Polygon using ArcGIS Desktop?

texts:
I have a number of polygons on which I wish to make analysis using points data symbols.

Is there a method through which I may add points automatically within each polygon in separate locations?

For Example I found centroid (x,y), and used it as one location. I need 12 locations in one polygon where the polygons are not of regular shape and points may not overlap too.


texts: ArcGIS JavaScript Extension for the Google Maps API, serverapi.arcgisonline.com/jsapi/gmaps/

texts:
My code using as the following:

Google Maps API in V3:

<script type="text/javascript" src="https://maps.googleapis.com/maps/api/js?libraries=geometry&v3&key=AIzaSyDYSz77yXE6n8okoIz4S16GbWvgPEyfUwk"> </script>

ArcGIS JavaScript Extension in V1.4:

<script type="text/javascript" src="http://serverapi.arcgisonline.com/jsapi/gmaps/?v=1.4"></script>

Then in my runtask(poly) function (which is adding polygon on the map):

mapExtension = new esri.arcgis.gmaps.MapExtension(map);
esri.arcgis.gmaps.Config.proxyUrl = "proxy.ashx";
var qtask = new esri.arcgis.gmaps.QueryTask("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/0");

But, my web application prompts an error like below:

  ?v=1.4:2 Uncaught ReferenceError: GLatLng is not defined_57 @
  ?v=1.4:2_256.execute @ ?v=1.4:2runtask @
  Default.aspx?facilityname=Airgas+National+Welders+-+10+&street=1221+New+Savannah+Road+&city=Augusta…:635mapplume
  @
  Default.aspx?facilityname=Airgas+National+Welders+-+10+&street=1221+New+Savannah+Road+&city=Augusta…:599initialize
  @
  Default.aspx?facilityname=Airgas+National+Welders+-+10+&street=1221+New+Savannah+Road+&city=Augusta…:452onload
  @
  Default.aspx?facilityname=Airgas+National+Welders+-+10+&street=1221+New+Savannah+Road+&city=Augusta…:991

I checked http://serverapi.arcgisonline.com/jsapi/gmaps/?v=1.4, it has GLatLng class which is under Google Maps API V2, but I am using V3 now, do you have new version working with V3 (google.maps.LatLng) ? I cannot find any versions applied to V3, but I want to keep using the same code, like esri.arcgis.gmaps.MapExtension function.


texts: Rasterio : rasterio.drivers() - 'module' object is not callable

texts:
I'm trying to pass a script on jupyter-hub with a specific kernel. With my last local installation it worked but now with my kernel I have this error:

TypeError                                 Traceback (most recent call last)
<ipython-input-27-f1a418595338> in <module>()
      4
      5 mask = None
----> 6 with rasterio.drivers():
      7     with rasterio.open(B2_reclass_path) as src:
      8         image = src.read(1) # first band

TypeError: 'module' object is not callable

I don't really understand this error and I didn't find a solution anywhere


texts: Exporting issue with colours in AutoCAD Civil 3D

texts:
I hope someone can tell me what's going on here.  On the left is what the image looks like in CAD and on the right is the result when I export to pdf format.  I also get the same color results when I send it to the printer.  As you can see, none of my other hatch patterns are giving me this issue besides the black area under "YOU ARE HERE".

Suggestions?  Thanks in advance.

ADDITION TO ORIGINAL QUESTION  I managed to achieve the appearance of the image on the left by adjusting the RGB value of the layer (very slightly).  The .pdf export and print both worked.  I changed one of the numeric settings such that the colour was not set to "absolute" black and it worked.  I would still like to know the answer to the original question though.  Although I've found a work-around, I still want to know why black does not print accordingly.


texts: Difference between WGS_1984_UTM_33N and GCS_WGS_84

texts:
What is the difference between:

a map that has a projected coordinate system of WGS_1984_UTM_33N, a
geographic coordinate system of GCS_WGS_1984, and a datum of
D_WGS_1984
a map that has no projected coordinate system, but with the same geographic
coordinate system of GCS_WGS_1984, and the same datum of D_WGS_1984


texts: How to add an existing feature class to a new feature dataset

texts:
I have a feature class.

IFeatureClass AOIfc = AOIfws.OpenFeatureClass("BrigadeAOI");

I need to be able to put this feature class into a feature dataset that is newly created.  The following line of code doesn't work for me.

IFeatureDataset AOIfd = AOIfc.FeatureDataset;

How do I accomplish this?  Thanks!


texts: Opening/saving QGIS project .qgs in different version

texts:
I was working on the latest QGIS 2.10 Pisa version and saved my work in .qgs project, I tried to re-open the projects with QGIS 2.0.1 Dufour but I had a lot of visualization problem and everything was messed up!

Any hints about it?

Iis there a way to save or convert a qgs project in older versions?


texts: Calculating distance from polygon to point

texts:
t1 = Polygon([(4.338074,50.848677), (4.344961,50.833264), (4.366227,50.840809), (4.367945,50.852455), (4.346693,50.858306)])
t = geopandas.GeoSeries(t1)
t.crs = {'init':'espg:4326'}

t2 = geopandas.GeoSeries([Point(4.382617,50.811948)])
t2.crs = {'init':'espg:4326'}

How can I calculate the distance of a point from the polygon respectively find the closest point between the polygon and the point?


texts: Plot two layers (raster and polygon-vector) with ggplot and use the fill-scale* for both layers

texts:
I am a little lost. I have a stars object where I would like to use the first attribute and map it to the
scale_fill*-function in ggplot.
I also have an polygon-object of class sf that has an attribute that I would also like to map to the
scale_fill*-function.
I really tried to compile some data and a reproducible example and I got this far:
It does not make any sense at all, but serves for the purpose:
library(stars)
library(tidyverse)
library(raster)
library(osmdata)
library(scico)

# get elevation
alt <- getData('alt', country='CHE')
alt_stars = st_as_stars(alt)

ch_canvas = st_bbox(alt_stars) %>% opq()

# overpass query
ch_moto = ch_canvas %>%
  add_osm_feature(key = "highway",
                  value="motorway") %>%
  osmdata_sf()

# get the sf object
df = ch_moto %>% pluck("osm_lines") %>%
  dplyr::select("maxspeed") %>%
  filter(!is.na(maxspeed),
         maxspeed != "none") %>%
  mutate(maxspeed = as.numeric(maxspeed))

# make the polygons
df_buffer = st_buffer(df, 0.1)

ggplot() +
  geom_stars(
    data = alt_stars,
    aes(x=x, y=y, fill = CHE_msk_alt),
    downsample = 2
  ) +
  scale_fill_continuous(low="white", high="red") +
  geom_sf(
    data = df,
    aes(fill=maxspeed)
  ) +
  scale_fill_scico(palette="berlin")

I am then getting this message:
Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.
I would just like to know how to plot any two layers (in many cases a stars-object and a sf-object) on top of each other and fill both according to the values in a variable.
I tried to use the cowplot-package and just define two ggplot objects and plot them above each other, but did not manage to make it work.
Maybe someone has a pointer on what to do in these situations:)


texts: Efficient spatial joining for large dataset in R

texts:
I am working with rather large data frames that I often need to do a spatial join on. The fastest way I have come up with so far is this method:

library(rgdal)

download.file("http://gis.ices.dk/shapefiles/ICES_ecoregions.zip",
              destfile = "ICES_ecoregions.zip")
unzip("ICES_ecoregions.zip")

# read eco region shapefiles
ices_eco <- rgdal::readOGR(".", "ICES_ecoregions_20150113_no_land", verbose = FALSE)
## Make a large data.frame (361,722 rows) with positions in the North Sea:
lon <- seq(-18.025, 32.025, by=0.05)
lat <- seq(48.025, 66.025, by=0.05)
c <- expand.grid(lon=lon, lat=lat)

# Get the Ecoregion for each position
pings  <- SpatialPoints(c[c('lon','lat')],proj4string=ices_eco@proj4string)
c$area   <- over(pings,ices_eco)$Ecoregion

But this takes a very long time and uses a lot of RAM, and will sometime come up with the Error: cannot allocate vector of size 460 Kb (if you can't reproduce the error, just make c larger).

Anyone can come up with a better/faster/more efficient solution?


texts:
I've created a bunch of Placemarks and have them in a kml file.
Is there a way to get this into the GPS of a car?


texts: QGIS Desktop 2.2.0 - program does not start

texts:
I just installed QGIS 2.2 (I download QGIS Standalone InstalleVersion 2.2 (32 bit)).Installation was successful without any errors.
Then I try lunch QGIS Desktop 2.2.0 but the program does not start.
Update
In postinstall.log I not found any bugs.
I also install  Microsoft Visual C++ 2008 Redistributable Package (x86) and clear path (Why is QGIS not launching?)


texts: Sample code to plan a trip using OpenTripPlanner as a jar dependency

texts:
I would like to use OpenTripPlanner as a jar dependency in my project and do some routing with it. I have followed the steps at https://github.com/opentripplanner/OpenTripPlanner/blob/master/src/test/java/org/opentripplanner/GtfsTest.java and they do seem to work if the start and end points of a RoutingRequest are near a transit stop but don't work otherwise.
 The following line returns me an empty list:

List<GraphPath> paths = new GraphPathFinder(router).getPaths(routingRequest);

I think its because the road network is missing from the graph. Specifically, the code in https://github.com/opentripplanner/OpenTripPlanner/blob/master/src/test/java/org/opentripplanner/GtfsTest.java only adds a transit network to the graph (the GtfsBundle). The road network is never added to the graph. So my question is could anyone please show me how to add the road network to the graph? I have the road network as an OSM PBF file. I seek a complete end-to-end example.

This link is saying at one point:

  once the data have been loaded the geometries of the  streets are
  extra information purely for presentation purposes.

Is the OSM road network not used for routing?


texts: Determination of point id after sampling in Google Earth Engine

texts:
Let's say I have a landcover map, assigning landcover codes to each pixel in a scene. Then I have a polygon and want to use some or all of the points in it to pick the landcover from the above map, and then sample a stack of Landsat images. I want to have control over final feature id's for my processing but I can not figure out how the id's are generated in GEE.

In the first stage, I will use either sampleRegions function (to retrieve all points in the intended polygon with their landcovers) or stratifiedSample function (to sample a fixed number of points for each landcover from the base map in the intended polygon). Base map is named 'image' and intended polygon is named 'geometry' i below code.

var points1 = image.sampleRegions({
          collection: geometry,
          scale: 30,
          geometries: true})
var points2 = image.stratifiedSample({
          numPoints: NumPoints,
          region: geometry,
          classBand: 'b1',
          geometries: true
          })

At this point, I get my point features with id's being integer numbers as 0, 1, 2, ... in stratifiedSample case, or the integers prefixed with with 0 (0_0, 0_1, 0_2, ...). What does this added '0_' mean? It is the same for all points, even if I have multiple polygons passing to sampleRegions.

Then I build a Landsat imageCollection for my time period of interest, and map it with sampleRegions function by passing any of above point sets to it. This time, a suffix (not prefix) of 0 is added to the point id's. So it will be either 0_0, 1_0, 2_0, ... (for stratifiedSample output) or 0_0_0, 0_1_0, 0_2_0, ... (for sampleRegions output).

var collection = LandsatCol
  .filterDate(start, end)
  .filterBounds(points)

var series =
    collection.map(function(image) {
      return image.sampleRegions({
        collection: points,
        properties: ['b1'],
        scale: 30,
        geometries: true
      });
    }).flatten()

I am confused what is the rule.


texts: How to set up future-proof line labels for data that is regularly updated?

texts:
I have a layer that I have processed heavily in FME before using it in QGIS. In the future I will have to reprocess it periodically when the source data is updated. It's a municipal street layer, which gets updated regularly.

The problem is that I require every street to be labelled in a destined-for-print map, which means a lot of manual editing of the labels to make it pretty. Since the label changes are stored within the table they will be overwritten when they table gets updated.

I'm looking for suggestions on how to avoid the problem.

I can think of 3 right off the bat, but don't know how to implement..

Don't store the label changes in the table. I have no clue how to approach this, but maybe it's possible in QGIS to use a table join and save the x,y, rotation channels in an external file?
With enough effort you can accomplish anything in FME, so presumably it's possible to update the table rather than overwrite it. My workspace is already fairly complicated though, and includes 2 linejoiners and 3 python scripts, so I can't imagine (yet) how to implement that.
Use a more sophisticated labelling tool that won't require so many edits. Labelling for print is a weakness in lots of GIS tools.

Option 1 seems simpler, but may not be possible.


texts: Circles and Projection in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2 and when I run topology on my feature class (which is in NAD 83 Lambert_Conformal_Conic) there were no overlaps. Once I project it into any projection I get overlapping polygons only where there are circles drawn with the circle tool.

This data will be provided to many customers that may choose any projection. The only way I have found to solve the problem is to go back and manually draw the circle instead of using the circle polygon tool. I know some quick tools to get rid of sliver overlaps, but I can't give this data to customers when they may choose many different types of projections. Is there another method to fix this problem?


texts:
The PostGIS documentation states that ST_PointOnSurface returns "a POINT guaranteed to lie on the surface". It seems like this function could be trivially implemented to give results that satisfy the documentation but provide little real-world utility, though I'm certain that PostGIS provides a non-trivial implementation.

This introduction to PostGIS provides a nice comparison and contrast of ST_Centroid with ST_PointOnSurface and says that "[ST_PointOnSurface] is substantially more computationally expensive than the centroid operation".

Is there a more thorough explanation of how ST_PointOnSurface is calculated? I've been using ST_Centroid, but have encountered some edge cases in my data where the centroid is outside the geometry. I believe that ST_PointOnSurface is the correct substitute, but the function name and documentation leave room for uncertainty.

Further, is the computational expense of ST_PointOnSurface incurred even if the centroid does lie within the geometry already?


texts: Why QGIS never uses more than 4GB of RAM

texts:
I see QGIS never uses more than about 4GB of RAM. I'm aware of the 32/64 bit difference, and in order to exclude some issues with my PC I tried a rather big Windows Server 2019 instance on AWS with 64GB of memory.
I had a very big raster file and I made some operations on it, but again I never saw RAM usage go above 4GB, the CPU itself seemed stuck at only 25%.
Is there anything that someone should tweak to unblock all the system resources?


texts: How do I show a popup on hover in ESRI Maps?

texts:
How do I show a popup on hover of a marker symbol in ESRI Maps?

Please note that I am not using leaflet JS.
like in the picture. its a demo from
https://embed.plnkr.co/393lgE49Ndqsj7O6dg19/
I need to find out whats equivalent of on('click/hover : ') for ESRI maps.


texts: Adding condition to relationship in ArcGIS Desktop?

texts:
I have the following data model:

Each point and each line can be linked to many measurements (1-M relationship)

I would like to add a condition in the relationship between MY_LINE and MEASUREMENTS, the relationship should only be active if A_BOOLEAN is TRUE

Is the best way to do that:

Add a new ID to MY_LINE and MEASUREMENTS which concatenate the fields THE_ID and A_BOOLEAN and then rebuild the relationship on this new ID;
Add a condition in the relationship somehow; or
something else?


texts: Adding the map canvas center to a composer label in QGIS

texts:
I am able to access the map canvas center in the python console in QGIS using

qgis.utils.iface.mapCanvas().center()

What I would like to do is get this into a map composer label.

Is this possible?


texts: Using Terrain Profile plugin with an interpolated raster

texts:
I am very new to QGIS and the goal of this exercise is to view an elevation profile of a river from a set of discrete points taken from surveying records.

Steps taken:

A vector layer was created from a delimited text layer (exported
from excel)
An .asc file was created using the interpolation plugin
The .asc file was converted to a raster .tif using the translate
        tool
Attempted to use the terrain profile plugin

I tried using a temporary polyline, but each time the plot showed up empty

Updates:

I got the following Python error when I tried adding my raster elevation layer to the Terrain Profile tool

@Bernd V. - I tried using a normal polyline, and still now success - see image below:


texts: QGIS placing points perpendicular on the travelling direction

texts:
Using QGIS 2.16.1 I have a line and a point shapefile. The line was made by using "points to route". The points have timestamps, in this case showing the direction of travel. I want to create a new point layer (not just labels...) perpendicular on the direction of travel. The practical point is to record the side of placing waste bins. If you are interested in the case: http://maasmission.wordpress.com

The points should only be placed on the right side (with an offset) of the direction of travel.


texts: Preventing clipped labels from QGIS Server WMS in Leaflet JS

texts:
I'm trying to show in Leaflet a WMS map served with QGIS Server. As you can see in this fiddle below there's a problem with the items in the map; they appear to be clipped somehow. The data which are given by QGIS Server are geographical multipolygons with associated their id in the postgis database.

https://jsfiddle.net/maurizioricci95/vydseLzq/1/

The map is working correctly if imported in QGIS Desktop using directly the WMS URL (http://quegis.labcd.unipi.it/cgi-bin/qgis_mapserv.fcgi). I also tried in Leaflet to request image in PNG format with the transparency but didn't fix the problem.

This is the GetCapabilities of the WMS:

http://quegis.labcd.unipi.it/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMTS&VERSION=1.3.0&REQUEST=GetCapabilities , everything appears to be ok.

It looks like there's some strange behavior/options set in Leaflet but I am not able to figure it out.


texts: Where can I get an OpenStreetMap extract that includes Bike Lanes?

texts: I'd like some way to download an extract of OSM data, like GeoFabrik or BBBike, that includes bike lanes. Does anyone have any advice or should I just try and figure them out with Osmosis?


texts: Resampling air temperature gridded data to sub-grid resolution

texts:
I'd like to obtain a 30-m resolution gridded dataset of air temperature from the ~1km gridded data of WorldClim, of Hijmans et al., 2005) for a mountainous region.

The easiest procedure (I thought it was...), with any Raster calculator, is T30m = T1km + (Z30m - Z1km)*lapse.rate, but this produces a 30-m grid on which the footprint of the 1km grid is kind of superimposed (map here: the colored cells "seem" to be at 1km but they are actually at 30m resolution, and you can see some small unexplained holes in the temperature raster)...

Does anyone have an idea to solve the problem, preferably with QGIS and/or R? Here is a sample of  the input data I use (30-m ASTER DEM, 1km GLOBE DEM and 1km air temperature).


texts: Displaying a box with text on a leaflet map

texts:
Is there a way to display a label or textbox on a specific geo-place on a map produced by leaflet.js.?

I know about the Marker- and Polygon Method But I can't find a way to put text in the polygon. I want the text and box to be the same size on every zoom level. And when dragging the map the box should be stuck to the map; like a sticky note.

Is there a way to achieve this with Leafletjs or with SVG?


texts: Speed up point sampling with ST_Value function in PostGIS, raster/vector overlay

texts:
I managed to extract raster values in PostGIS based on a table of point geometries by following Point sampling raster with PostGIS

SELECT id, ST_Value(rastertable.rast, 1, pointtable.geom) as RasterValue
FROM rastertable, pointtable
WHERE ST_Intersects(rastertable.rast, pointtable.geom)
LIMIT 20

The above query takes roughly 20 seconds, which is one second per case. In the end I want to extract raster values for over 10 million points. With the current query that would take roughly 115 days; not really an option. How could I adjust the above query to increase the performance? Or do I maybe even need to change the method?

Edit:
Using QGIS point sampling plugin also results in a 3 second per 10 point speed, reducing the waiting time to about 40 days. That still is not workable.

Version info: