texts: How to get counts from a table to produce a heatmap?

texts:
I am trying to create a heat map of doctors per zip code based on zip code. I have a table with the doctors and their zip codes, some of them repeating because there is more than one doctor per zip code. I also have a shapefile with polygon features of zip codes.

I have tried the Get Count tool (which returned the total number of rows) and the Find Duplicate tool, which did nothing for me. Find Duplicate just produced a confusing table with the same amount of records as the original that did not actually have the correct input values in the center column (many were something different) and the counts in the right column did not match the counts I see (via select by attribute: zipcode = thiszipcode, result: x were selected).

My thought is that I could create a table based on the doctors table that lists each unique value (the zipcode) and its count. I would then join that table to the zip code shapefile and create a nice heat map. It seems so simple, but I cannot figure out how to get the count. Do you have any ideas?


texts: Using a Variable in arcpy.SelectLayerByAttribute

texts:
I know this is a very basic problem I should be able to solve myself, but I've had no luck. I have a variable I want to use in the query for arcpy.SelectLayerByAttribute, but the line that was working yesterday is now giving me an "error 000358: invalid expression". I simply want to use the variable CountyNM in the where clause.

Here is the section of code giving me issues:

CountyNM = "Hamilton"
arcpy.SelectLayerByAttribute_management("hzCountyLYR", "NEW_SELECTION", 'CountyName =' + "'"+ CountyNM +"'")

Syntax: SelectLayerByAttribute_management (in_layer_or_view, {selection_type}, {where_clause}, {invert_where_clause})


texts: Color feature according to datetime field

texts:
I have a layer in QGIS which comes from a CSV with almost a million points. One of the columns of this csv is a datetime.

I would like the color of the points to be on a gradient depending on this datetime (e.g. oldest points in blue, newest in red, continuous gradient). I'm new to QGIS but know about the "Categorized" marker type, which allows me to have different markers depending on a field. Unfortunately I can only choose numerical fields there, and not text or datetime fields.

Is there a way to do this?

Alternatively, I could sort the source CSV by date. Is there a way to then use the "line number" to do the color gradient?

My CSV data looks like this:

Accuracy,Altitude,Bearing,Date Time,Date(GMT),Description,Distance,Latitude,Longitude,Name,Speed,_FileNames,Location,NAME_0,NAME_1,NAME_2,NAME_3,NAME_4,NAME_5,REGION,UID,Time elapsed
13.0,3089.800049,12.0,2013-03-16 18:03:00,,,24156.0,-13.322358,-72.049927,,10.5268,C:/Users/Alexis/Dropbox/Location/Aggregated.csv,Null,Peru,Cusco,Urubamba,Huayllabamba,,,,135966.0,0.0
26.0,2941.199951,17.20000076,2013-03-16 18:03:00,,,24156.0,-13.321477,-72.049874,,2.5,C:/Users/Alexis/Dropbox/Location/Aggregated.csv,Null,Peru,Cusco,Urubamba,Huayllabamba,,,,135966.0,0.0

Note that I have changed in the layer's properties-->fields, I have changed the field type to datetime, according to the datetime format in my data.


texts:
I have trouble categorize my markers by category (boolean in this case) in geoserver using CSS.
I have tried a few different combinations without success.

Simply I want: if true - yellow square, if false - blue

Suggestions?

What I have tried:

* {
  mark: symbol('square'), symbol('square');
}

    [exist=true] :nth-mark(1) {
      size: 12px;
      fill: yellow;
      stroke: yellow;
    }

    [exist= false] :nth-mark(2) {
      size: 12px;
      fill: blue;
      stroke: blue;
    }

[exist=true] {
  mark: symbol('square');
  size: 12px;
  fill: yellow;
  stroke: yellow;
}

[exist= false] {
  mark: symbol('square');
  size: 12px;
  fill: blue;
  stroke: blue;
}

If I may ask in the same question.
This is from doc.geoserver, what is usage, column name?

[usage='hospital'] {
  mark: symbol('circle'), symbol('cross');
}


texts: How to read csv with WKT geometry in Arcmap?


texts: ArcGIS Server 10.1 GP Service non-functional after ArcSDE restart

texts:
every night, the Oracle/ArcSDE server gets restarted. As a result, some, but not all, map and GP services require a restart to be functional afterwards. Can others confirm this experience and share solutions?

at 10.0, when publishing a service, there was an option for the service to check for this sort of thing and repair itself, but this seems to have disappeared at 10.1.

failing and other feedback, we propose to automate the restart of all ArcGIS Server services after the nightly ArcSDE restart.


texts: arcpy 'Feature Layer' conversion to 'KML file' doesn't do anything (doesn't even throw error)

texts:
I have built a small arcpy script to convert feature classes into kml files (kmz in fact).

import arcpy
import os
import datetime

working_folder = os.path.normpath(os.path.join(os.path.dirname(__file__), '..'))

#Creating the path where the output kmz files will be written
date = datetime.date.today()
out_folder_name = 'Wibre Network ' + str(date)
out_path = os.path.join(working_folder, 'DATA\KML files', out_folder_name)

#Input variables.
fc_Opt_All = os.path.join(working_folder, 'DATA\Optimity_Equipment.gdb', 'Optimity_All')

#feature class path
fc_node_live_available = os.path.join(working_folder, 'DATA\Optimity_Equipment.gdb', 'Node_Live_Available')

#kmz file path
kml_node_live_available = os.path.join(out_path, 'Node_Live_Available.kmz')

#Conversion
arcpy.MakeFeatureLayer_management(fc_Opt_All, 'fl_Opt_All')
query = "ESTATE LIKE '{0}' AND STATUS LIKE '{1}'".format('Node Live','YES')
arcpy.SelectLayerByAttribute_management('fl_Opt_All', 'NEW_SELECTION', query)
arcpy.CopyFeatures_management('fl_Opt_All', fc_node_live_available)
arcpy.MakeFeatureLayer_management(fc_node_live_available, 'fl_node_live_available')
arcpy.LayerToKML_conversion('fl_node_live_available', kml_node_live_available, is_composite = 'NO_COMPOSITE', ignore_zvalue = 'ABSOLUTE')

The script doesn't throw any error. However, it doesn't produce any output.
I have double checked the feature class created after the query (fc_node_live_available) and it's NOT empty.

Any ideas of what's happening?


texts:
Hello I am new to working with GeoTiff data and I am trying to grab lat/lon values from a tiff image. Here is my code:

import gdal
import matplotlib.pyplot as plt
import osr

# open the dataset and get the geo transform matrix
ds = gdal.Open(r'C:\Users\_______\Documents\Landsat8\LC08_L1GT_156120_20170207_20170216_01_T2_B5.tif')
xoffset, px_w, rot1, yoffset, rot2, px_h = ds.GetGeoTransform()

x = 100
y = 100

# supposing x and y are your pixel coordinate this
# is how to get the coordinate in space.
posX = px_w * x + rot1 * y + xoffset
posY = rot2 * x + px_h * y + yoffset

# shift to the center of the pixel
posX += px_w / 2.0
posY += px_h / 2.0

# get CRS from dataset
crs = osr.SpatialReference()
crs.ImportFromWkt(ds.GetProjectionRef())
# create lat/long crs with WGS84 datum
crsGeo = osr.SpatialReference()
crsGeo.ImportFromEPSG(4326) # 4326 is the EPSG id of lat/long crs
t = osr.CoordinateTransformation(crs, crsGeo)
(lat, long, z) = t.TransformPoint(posX, posY)

However, when I run this I keep receiving the following error:

NotImplementedError: Wrong number or type of arguments for overloaded
function 'CoordinateTransformation_TransformPoint'.
Possible C/C++ prototypes are:
OSRCoordinateTransformationShadow::TransformPoint(double [3])
OSRCoordinateTransformationShadow::TransformPoint(double [4])
OSRCoordinateTransformationShadow::TransformPoint(double [3],double,double,double)
OSRCoordinateTransformationShadow::TransformPoint(double [4],double,double,double,double)

I have no idea what this means, I have looked up the error and found no useful information so I am hoping one of you have experienced this problem before.


texts: Using difference/clip tools in QGIS without it creating blank layer?

texts: I am using QGIS in an effort to find addresses without municipal water in arkansas.  I have a polygon layer of Arkansas public water systems as well as a point layer of arkansas residencies.  I believe I need to use the difference tool to be able to cut out all of the residencies that are overlapping with the public water systems, but everytime I try that, it just creates a blank shapefile layer. How should I go about cutting the layers?


texts: How to vectorize a Raster-Classification? (QGIS, GRASS)

texts:
I classified my raster-DHM into elevation levels using colormap. Its not a "real" classification, just a colormap. So I can't vectorize the classes to get polygons for every class.

What's the best way, to classify and vectorize a raster? The Rastercalculator? How to use it, if I want to classify my DHM in elevation gradation like: 1000-1200m; 1200-1400m?


texts: Need second-level subdivisions (Counties) for the United Kingdom

texts:
I am working on a map that requires second-level subdivisions for the United Kingdom, looking for these in shapefile format. These would be county-level subdivisions.  I've done a lot of searching at the data.gov.uk website and also at Ordnance Survey GIS data download website but I can' t seem to find the right shapefile.

Ideally the shapefile should have the second level subdivisions that match the current Wikipedia resource at
https://en.wikipedia.org/wiki/ISO_3166-2:GB#Second-level_subdivisions,
so  Buckinghamshire should have the code GB-BKM in the shapefile, Cambridgshire should have the coding GB-CAM   Cambridgeshire, and Cumbria should have the coding GB-CMA.

Does anyone know of a download link from data.gov.uk or another UK government website that has this? I've done numerous online and Google searches. I found this one one on ceremonial counties - https://data.gov.uk/dataset/ceremonial-county-boundaries-of-england but it doesn't seem like the right one.  And I'm Canadian, so I'm confused by the different offerings out there such as Ceremonial counties, and historical counties, unary, etc.


texts: How to implement 4D Choropleth Map using cesiumjs?

texts:
I want to implement 4D Choropleth Map on cesiumjs. Anyone know how to do it and show me the code?

(source: cesiumjs.org)

(source: cesiumjs.org)

(source: cesiumjs.org)


texts: Normalize large rasters with gdal.translate

texts:
i have a series of ndvi images that i am trying to normalize with gdal to run some calculations, but i think i'm missing something.

per:
    https://github.com/OSGeo/gdal/blob/master/gdal/swig/include/python/gdal_python.i

scaleParams --- list of scale parameters, each of the form [src_min,src_max] or [src_min,src_max,dst_min,dst_max]

my original  stats are

[ STATS ] =  Minimum=-0.194, Maximum=0.459, Mean=0.103, StdDev=0.073

i set my newmin and newmax to 0 and 1

then i

my_options = gdal.TranslateOptions(scaleParams = [[-0.194, 0.459,newmin, newmax]], format = 'GTiff',
                                      creationOptions = ['TFW=YES', 'COMPRESS=LZW']

gdal.Translate(gdaloutput, my_raster, options=my_options)

yet my output is

[ STATS ] =  Minimum=-30.148, Maximum=105.483, Mean=26.477, StdDev=18.837

I am getting the statistics with

srcband = gtif.GetRasterBand(1)

# Get raster statistics
stats = srcband.GetStatistics(True, True)

I had thought this was the correct way to get the statistics?
https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#loop-through-all-raster-bands

i tried with a different raster and changed my min and max to floats, from 0 and 1 to 0.0 and 1.0. which gave me

[ STATS ] =  Minimum=0.226, Maximum=1.083

from the original of

[ STATS ] =  Minimum=-0.211, Maximum=0.439

I also read somewhere about avoid setting the new max higher or the new min lower than the original, so following that, setting new min to 0 and new max to 0.439, for the new raster

i got

[ STATS ] =  Minimum=0.099, Maximum=0.475

which places my new max higher than the existing.

What am i doing wrong?

Edit:

without the double brackets

Changing

[[src_min,src_max,dst_min,dst_max]] to [src_min,src_max,dst_min,dst_max]

you get

TypeError: 'float' object is not iterable

both

ComputeRasterMinMax()  and ComputeStatistics()

returns

(-0.19442784786224365, 0.4585714042186737)


texts: Finding administrative boundaries in Guinea using OSM

texts:
I am trying to map the quarters (neighborhoods) in the city of Conakry, Guinea.

As far as I can find, OMS contains information on the location of each quarter, and displays this with a point in the center of each quarter.

I need the boundaries of each quarter. I have checked www.gadm.org, they have boundaries up to the point of commune (sub-prefecture), which is one level above quarter.

My question is, is this information in OSM, and did I miss it (no purple lines), if not, is there an other way to find it?


texts: How to export color of polygon from Autocad into QGIS?

texts: I use Qgis 2.14. I have some problem with export .dgn to shp from Autocad. Qgis appropriates random colors to my polygons after export. Somebody knows how i can export original color of my Autocad dgn?


texts: Trying to implement Albers projection in Python, what is going wrong?

texts:
I have a platform where I can't import Proj4. I only need to convert points to Albers projection, so it seemed simple, but my conversion has some errors:

Projection definition: +proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105 +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs'

Here's the code implemented from the formulae on Wikipedia:

from math import sin, radians, cos

def albers_proj(lon0, fi0, fi1, fi2, R=6367444.5, x0=0):
    lon0, fi0, fi1, fi2 = [radians(i) for i in (lon0, fi0, fi1, fi2)]

    n = (sin(fi1) + sin(fi2)) / 2
    C = cos(fi1) ** 2 + 2 * n * sin(fi1)
    ro0 = R / n * (C - 2 * n * sin(fi0)) ** .5
    def transformer(lon, lat):
        lon, lat = radians(lon), radians(lat)
        ro = R / n * (C - 2 * n * sin(lat)) ** .5
        theta = n * (lon - lon0)

        return ro * sin(theta) + x0, ro0 - ro * cos(theta)

    return transformer

albers_siberia = albers_proj(105, 0, 52, 64, R=6367444.5, x0=18500000)
print(albers_siberia(83, 55))

Output:

(17127237.118954152, 5777762.666018201)

I found the best R to minimise the error, but still something is not right:

    wgs84: 83.0 55.0

r = 6356302.97445488

     qgis: 17121825.913615   5767652.8720776
my script: 17129639.13010738 5767652.944553478
    error:     7813.21649238       0.072475878

What am I doing wrong?


texts: Disable 'View Feature Details' panel on geocortex

texts: Does anyone know if its possible to disable the ability to drill down from the initial results panel to the additional feature details panel? My feature description has all the feature attributes I'd like to display so wish to disable the next step. I have made the attributes invisible so the Additional Feature Details are blank. .


texts: Links will not open to new window in Carto?

texts:
I would like to have my links open in new windows in CartoDB which I thought I had achieved successfully in the past using the following HTML in the infowindow.

STREET VIEW

This code no longer seems to work on my new maps and even my old maps. Has anyone else encountered this? What HTML do you use to open in new tabs, windows, etc.

I've noticed that when just adding the url field by itself with no HTML formatting it will work in CartoDB editor (even opening in a new window) but is no longer recognized as a URL once embedded.

Below is a sample google street view I would be linking to:

https://www.google.com/maps?q&layer=c&cbll=40.7307132,-75.313928,&cbp=11,0,0,0,0


texts: Calculate the line length between consecutive points in a network in arcgis

texts:
I am trying to find the road distance between consecutive or say all nearest points.

A sample of my network is shown in the image. My actual network has more than 5000 points. Now I want to calculate the distance between consecutive points like point id 335 and 3797, then 3797 and 2281 and so on, the disconnected points can be ignored.

I tried using closest facility option. I loaded all the points as facilities and incidents to find 3 nearest facilities but it is giving me around 40k + routes which are repeated/overlapping and also the longer routes are not covered, like the nearest facility from point id 2281 will be 2281, 3797 and 678 and it will not cover 1818

I am not sure if i have started in the right direction. If anyone has solved similar kind of problem kindly help me out.

Notes about the network: all the points are properly placed, there are no overlapping lines, no snapping problem, the topological errors have been removed.


texts: How to fix a "arguments did not match any overloaded call" error

texts:
I'm making a geological map with the use of SpatialLite points, lines and polygons, but suddenly QGIS showed all kind of errors and decided it doesn't want to load my maps and SpatialLite layers any more. It happened when I clicked on a map layer to turn it on.

This is what it says:

  TypeError: arguments did not match any overloaded call:   QSize(): too
  many arguments   QSize(int, int): argument 1 has unexpected type
  'long'   QSize(QSize): argument 1 has unexpected type 'long'
  Traceback (most recent call last):   File
  "C:/Users/kikid/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py",
  line 125, in request
      self.setup_map()   File "C:/Users/kikid/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py",
  line 201, in setup_map
      olSize = QSize(int(olWidth), int(olHeight)) TypeError: arguments did not match any overloaded call:   QSize(): too many arguments
  QSize(int, int): argument 1 has unexpected type 'long'   QSize(QSize):
  argument 1 has unexpected type 'long'

  Python versie: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64
  bit (AMD64)]  QGIS versie: 2.18.6 Las Palmas, 5f369b0

  Python pad:
  C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing
  C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python C:/Users/kikid/.qgis2/python
  C:/Users/kikid/.qgis2/python/plugins
  C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins
  C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
  C:\PROGRA~1\QGIS2~1.18\bin C:\PROGRA~1\QGIS2~1.18\apps\Python27
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\PIL
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
  C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
  C:/Users/kikid/.qgis2//python C:/Users/kikid/Desktop/Kiki BSc

I tried to load my SpatialLite layers in a different QGIS map and that worked, it only didn't contain all data and colours I gave it, so it will be a loooot of work to fix that all.

Is there any way to get rid of this error?


texts: Problem with insert into PostGIS view table via QGIS

texts:
I use a PostgreSQL 9.5.14 (PostGIS) and QGIS 2.18.7.

The problem is that when I try to add a new object to the View via a QGIS, the following error occurs:

  ERROR:  column "geom" is of type geometry but expression is of type integer
      LINE 1: INSERT INTO testtr VALUES (NEW.objectid, NEW.num_pdr, NEW.ge...
                                                       ^
      HINT:  You will need to rewrite or cast the expression.
      QUERY:  INSERT INTO testtr VALUES (NEW.objectid, NEW.num_pdr, NEW.geom)
      CONTEXT:  PL/pgSQL function testtr_view_dml() line 4 at SQL statement

But update and delete triggers work properly.
I think problem is in the trigger function for insert.

My query:

At first i create a table in pgAdmin with next query

-- Table: public.testtr

-- DROP TABLE public.testtr;

CREATE TABLE public.testtr
(
  id integer NOT NULL DEFAULT nextval('testtr_id_seq'::regclass),
  geom geometry(MultiPointZ,100000),
  num_pdr integer,
  CONSTRAINT testtr_pkey PRIMARY KEY (id)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE public.testtr
OWNER TO postgres;
GRANT ALL ON TABLE public.testtr TO postgres;

Next I create a view for this table with query

-- View: public.testtr_view

-- DROP VIEW public.testtr_view;

CREATE OR REPLACE VIEW public.testtr_view AS
 SELECT row_number() OVER (ORDER BY testtr.*) AS objectid,
    testtr.num_pdr,
    testtr.geom
FROM testtr;

ALTER TABLE public.testtr_view
OWNER TO postgres;
GRANT ALL ON TABLE public.testtr_view TO postgres;

-- Trigger: testtr_view_dml_trig on public.testtr_view

-- DROP TRIGGER testtr_view_dml_trig ON public.testtr_view;

CREATE TRIGGER testtr_view_dml_trig
  INSTEAD OF INSERT OR UPDATE OR DELETE
  ON public.testtr_view
  FOR EACH ROW
  EXECUTE PROCEDURE public.testtr_view_dml();

I need this view to be editable in QGIS so I create a trigger fuction

-- Function: public.testtr_view_dml()

-- DROP FUNCTION public.testtr_view_dml();

CREATE OR REPLACE FUNCTION public.testtr_view_dml()
  RETURNS trigger AS
$BODY$
    BEGIN
        IF TG_OP = 'INSERT' THEN
        INSERT INTO testtr VALUES (NEW.objectid, NEW.num_pdr, NEW.geom);
        RETURN NEW;

        ELSEIF TG_OP = 'UPDATE' THEN
        UPDATE testtr SET geom = NEW.geom WHERE id=OLD.objectid;
        RETURN NEW;

      END IF;
      RETURN NEW;
    END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION public.testtr_view_dml()
  OWNER TO postgres;

AND at last the trigger

-- Trigger: testtr_view_dml_trig on public.testtr_view

-- DROP TRIGGER testtr_view_dml_trig ON public.testtr_view;

CREATE TRIGGER testtr_view_dml_trig
  INSTEAD OF INSERT OR UPDATE OR DELETE
  ON public.testtr_view
  FOR EACH ROW
  EXECUTE PROCEDURE public.testtr_view_dml();


texts: How to curve lines between points on a map?

texts:
I'm trying to curve the lines between points on a RaphaelJS map like Tata Communications - but how did they do it? What sort of mathematics lie behind these curves - unless ofcourse they're all manually drawn?

Some dude told me Catmull-Rom Splines, but those will only smoothen lines, not curve them.

Many thanks!


texts: Select Layer By Location to feature class in ModelBuilder?

texts: What is the best way to convert the Feature Layer output of Select Layer By Location to Feature Class in ModelBuilder?


texts:
Hopefully someone can help. I'm using OpenLayers and have a TMS layer with local tiles stored. When drawing the tiles it seems slow, I see the tiles load one by one.

Is there a way to make getURL request handle multiple tiles at once?


texts: Error 000865 when using arcpy.ASCIIToRaster_conversion

texts:
I don't know what I'm doing wrong. I always get the error 000865 when using  arcpy.ASCIIToRaster_conversion. I'm new and do not have much Python experience right now – so sorry for that.

>>> import arcpy
import os
from arcpy import env
env.workspace = r"D:\Julia_T\projekt\erste_schritte\dritter_tag_python\testdaten"

#folder to list
ascFileList = arcpy.ListFiles("*.asc")

for ascFile in ascFileList:
# geoprocessing steps
ascFileName = os.path.splitext(ascFile)[0]
# define the output file
rastFile = ascFileName + "_output.img"
# run the tool
arcpy.ASCIIToRaster_conversion(ascFile,rastFile, "INTEGER")

  Runtime error  Traceback (most recent call last):   File "", line 15, in         File "c:\programme\arcgis\desktop10.2\arcpy\arcpy\conversion.py", line 2193, in    ASCIIToRaster     raise e ExecuteError: ERROR 000865: Eingabe-ASCII-Raster-Datei:    2001002_cloud_free.asc ist nicht vorhanden.

Do I have to use the CatalogPath? I'm working with ArcMap 10.2.1.
Thanks


texts: shapefile record has no attribute named bounds

texts:
This is the code that I am using to print the bounds of my shapefile.

import fiona
import shapely
import shapely.geometry

fc = fiona.open("filename.shp")
shapefile_record = fc.next()
shape = shapely.geometry.asShape(shapefile_record['geometry'])
polygon = shapefile_record['geometry']
print polygon.bounds

However, I keep getting an error that says:

AttributeError: 'dict' object has no attribute 'bounds'

Can someone guide me as to where I am wrong?


texts:
Does anybody know where I can find the technical details for Burma's latest datum? It's called Myanmar Datum 2000 and I can't find any details online. I need the details to create a custom datum on my Garmin GPS.

What I need is the dx, dy, dz, da, and df.

If this isn't available, what is the closest datum that could be used instead?


texts: Split Multiple Rasters into multiple rasters of different size in QGIS

texts:
I am wondering if there is a way to take multiple .tiff or .ecw files of any given size, and divide them up into multiple files of the same type, but of different size, using QGIS.

Also, is it is possible to do it by degrees, rather than "rows" and "columns"?

For example, say I have aerial imagery of a given region in the world. The total imagery is 10 degree by 10 degrees (lat/lon), but is divided into 4 equal quadrant images. Now say I want it to be divided into 1x1 degree (lat/lon) blocks. It would go from 4 image files (2x2) to 100 (10x10).

It seems like I could do this manually by first merging the rasters as explained here, then simply adding that single large raster to the canvas as selecting Save As and using the VRT Tile option as explained here to split it. I tried to use the merge capability to do this, but my PC locked up for hours. I couldn't move anything and ended up just hard powering it off. There was a very small corrupt file that I just trashed, so it didn't seem to get far.

My question is, do you even need to do the merge. Can you split multiple files? If not, I might just write a new plugin that will.


texts: Gdal to acquire Netcdf-like metadata structure in python

texts:
I'm working with gdal library after studying the Netcdf-4 one. Is there a way to produce a json metadata output similar to the one offered by the latter lib?

this is what i got so far:

from osgeo import gdal

gdal.UseExceptions()

file = gdal.Open("NETCDF:../files/ECMWF_ERA-40_subset.nc")

filemeta = file.GetMetadata()

subsmeta = []

for subs in  file.GetSubDatasets():

    subsmeta.append(gdal.Open(subs[0]).GetMetadata())

but the output is much more confusing than the relative "ncdump -h" version and doesn't include dimensions metas.


texts: Correcting forestry acreages due to slope to get actual acreages treated

texts: I came up with a process to get actual acreages vs map acres of steep slopes, but I am wondering what others do to correct for steep slopes and get actual acreage treated? For mine, In ArcGIS Desktop, I use a DEM to create a percent slope raster. After classifying the slope data, I use our treatment polygons to cut out the slope data so I know how many acres I have at a particular slope. Then I use a table found in a range management book that tells you what percentage to increase the acreage by per slope percent. It works well, I'm just wondering if there is a better more efficient process or tool I could be using?


texts: How to force a ol.source.Vector to reload data from server?

texts:
I have a ol.source.Vector like this:

var geoJSONFormat = new ol.format.GeoJSON();

vectorSource = new ol.source.Vector({
    strategy: ol.loadingstrategy.bbox,
    loader: function(extent, resolution, projection) {
        var url = getUrl( extent );
        $.ajax({
            url: url,
            dataType: 'json'
        }).then(function(response) {
            var features = geoJSONFormat.readFeatures( response, {featureProjection: 'EPSG:3857'} );
            console.log( features );
            vectorSource.addFeatures( features );
        })
    },
    format: geoJSONFormat
});

mapillaryLayer = new ol.layer.Vector({
    source: vectorSource
});

map.addLayer( mapillaryLayer );

When first load, all its ok. I can see the features on the map. When I pan the map to an undiscovered area, the features are updated from server.

But when I zoom in and pan to an already seen area (inside the original wide area before the zoom in) no new loads are made.

I know its a correct behaviour because all data inside this BBOX must be taken when zoomed out, so when zoom in I already have all I need for that area.

But I'm loading a limited amount of points in my query (maxresults):

var url = "getPhotosInBBOX?minlat="+minlat+"&minlon="+minlon+"&maxlat="+maxlat+"&maxlon="+maxlon+"&maxresults="+maxresults;

And now, even already loaded the BBOX content, I need to do a refresh because when zoom in, I need to have the same amount of data I have before.

Example: I have a bbox to show all Canada territory. So I will ask for 100 features from my server. Now, I'll zoom to Alberta. Alberta is inside Canada bbox and the layer will not load again. But I need to refresh to show 100 features inside this new bbox.

Why? My query is random. It takes 100 arbitrary features inside the view's current extent. May all be in Alberta, may not. When I zoom in, I need to take those I don't take before.

I don't know I'm clear. Ask for further informations.


texts: Extracting number of vertices in each polygon?

texts:
I have ArcGIS Desktop 10.2 and my challenge is how to extract the of number vertices in each polygon for all features like this:

I have many parcels in my feature class and I want to extract the number of vertices for all features separately, then I want to display XY coordinates for all vertices.

for more information, i just want to convert vertices for each polygon and display number of each vertices staring from 1 number, so if i have polygon and it have a 4 vertices, a want to convert polygon to vertices i display number of vertices like this (1,2,3,4,5), then displaying xy for each vertices, i think the real challenge id to, how to convert all polygon to vertices and make each number of vertices starting from 1 number.


texts: What are spatial and attribute indexes in geodatabases?

texts:
What is the difference between a spatial index and an attribute index in an ArcGIS geodatabase?

What are their purposes?


texts: Using NetCDF file in PostgreSQL database?

texts: I have a NetCDF file which is 25 by 20 grids by latitude and longitude and contains hourly rainfall data at each grid. In total I have 780 hours worth of data. I was trying to use ArcGIS to create a raster, but I could only get one hour at a time. Is there a way to convert it to another file type or to a raster with all the hours so I can easily use it in my database?


texts: Arcpy error message: 'Layer' object is not subscriptable?

texts:
What the code does successfully so far, is gets a list of the layers in the first data frame of the current map. What it doesn't succeed in doing is, joining the layer names to the gdb location using os.path.join(), so I can then do a clip against an existing feature class.

I think I need to convert the lyr names to a string before I can use it in a Clip_analysis, how do I do that?

Here is what I've got so far:

import arcpy
import arcpy.mapping
import os

#Set the current map
mxd = arcpy.mapping.MapDocument("CURRENT")

# Set the dataframe
df = arcpy.mapping.ListDataFrames(mxd)[0]

# Get list of Layers from toc
layers = arcpy.mapping.ListLayers(mxd,"",df)

#Set the output workspace
outWorkspace = r'C:\Users\jsommerville\Documents\ArcGIS\California\ProjectData.gdb'

for lyr in layers:
    outFeatureClass = os.path.join(outWorkspace, lyr)
    arcpy.Clip_analysis(lyr, "extent", outFeatureClass)

The error message I'm getting: 'Layer' object is not subscriptable


texts: GCS_Clarke_1880 transformation to GCS_WGS_1984 in ArcGIS 10.3.1

texts:
ArcGIS 10.3.1 has no transformation for GCS_Clarke_1880 to GCS_WGS_1984. Why is this the case, what have I missed?

Do we have to write a custom transformation?

There are a heap of Clarke transformation equations on stack. But I don't understand why.

GCS_Clarke_1880 WKID: 4034 Authority: EPSG

  Angular Unit: Degree (0.0174532925199433)  Prime Meridian: Greenwich
  (0.0)  Datum: D_Clarke_1880    Spheroid: Clarke_1880
       Semimajor Axis: 6378249.144808011
       Semiminor Axis: 6356514.966204133
       Inverse Flattening: 293.4663076556253

To

GCS_WGS_1984 WKID: 4326 Authority: EPSG

  Angular Unit: Degree (0.0174532925199433) Prime Meridian: Greenwich
  (0.0) Datum: D_WGS_1984   Spheroid: WGS_1984
      Semimajor Axis: 6378137.0
      Semiminor Axis: 6356752.314245179
      Inverse Flattening: 298.257223563


texts: ArcObjects ISQLSyntax with Oracle DB esriSQL_DelimitedIdentifierPrefix issue

texts:
New to stack exchange, heading over from the esri forums, it is great!

I am trying to produce a query that I can use across all workspaces, naturally using ISQLSyntax. The problem I'm having is that when using esriSQL_DelimitedIdentifierPrefix on a Oracle SDE instance it is returning a speech mark, whereas I'd expect it to return a blank string. When I use this code on a FileGDB or an MDB it works as expected, but on SDE I get the error:

Underlying DBMS error[ORA-00904: "table.field_name": invalid identifier][OWNER.TABLE]

The error goes away if I remove the prefix and suffix, so there are no typos in the query. I don't want to do a select case SDE because it seems against the point of having ISQLSyntax.

Any Ideas?

Dim Prefix As string, Suffix As string
Dim upper As String, sQualColName As String
Prefix = pSQLSyntax.GetSpecialCharacter(esriSQLSpecialCharacters.esriSQL_DelimitedIdentifierPrefix)
Suffix = pSQLSyntax.GetSpecialCharacter(esriSQLSpecialCharacters.esriSQL_DelimitedIdentifierSuffix)

sQualColName = pSQLSyntax.QualifyColumnName("table", "field_name")

sVal = Prefix & sQualColName & Suffix & " = " & Chr(39) & Trim(strFindString).ToUpper & Chr(39)


texts: Input data to Maxent niche modelling tool?

texts:
I have following question regarding the input format of files to Maxent.

From here I downloaded the environmental layers but they are in tiff format and I guess the tool needs in asc format where number of columns,rows etc and something like a matrix is there.  How do I get the exact format required by the tool?

How do I solve it.


texts: Is there an existing solution for pattern recognition?

texts:
I've been looking at modeling potential avalanche hazard based on terrain characteristics.  I have a vector layer showing historical avalanche paths and have raster layers for elevation, aspect, slope, curvature and land cover.

I was hoping there was an existing way to use the vector layer to identify what value range from the underlying rasters were present in the known avalanche paths and, perhaps using a bayesian statistical approach, predict/identify new areas that have similar underlying characteristics.

Does a validated GIS model or process already exist to approach this?


texts: how to add polylines(routes from public transport) on postgresql(pgrouting) with GUI

texts:
After searching a several long time... i decide to ask here for help,
  im currently using ubuntu 12.04(because of compatibility) with
  postgresql 9.1, postgis2, pgrouting2 and qgis2.2

i need to make a system for network adding one by one the lines and their "fixed routes" to the database, with eclipse/c# or something like that, so i can do custom interface to add polyliness into database with a GUI :

Polylines(routes with name tag or ID, for later use on queries)
A polyline is composed by 2 routes (one to go and on to back, where the last node of one is the first of the other like a circle)
Polylines have infinite reverse cost(one way direction).
Some linestrings might use some of the same nodes of other line.

I'm kinda newbie in this GIS world, so a detailed guide would be nice.

i read that each thread should contain exactly one question and its answers, so other parts would go in other threads.

EDIT2:

"lines", mean public transport routes
i have a network(streets) that i want to extend with public
transport (PT) routes in another layer
i want to add one by one because there is no PT info available(lines
or their routes)

EDIT3: i did some research again and see that what i need are polylines for each of my routes

i have a network(streets) like this:

i want to add polyliness like this:

and one by one like this (but in the same table/layer, so i can use them on shortest patch query later):


texts:
I have two tables in a personal geodatabase and one of them has 13 records and 6 fields, the other has 118 records and 125 fields.

There is a relationship class between the two tables where the first table is the origin and the other (bigger) table is the destination.

What I need is to take the records satisfying a condition from the second table and sum those records and put into a new field in the origin table


texts: Comparing various JavaScript mapping libraries?

texts:
I am working on a web-based mapping system and I'm trying to figure out which library to use.

These are links to comparisons of available libraries:

Laurent Jegou's benchmark (from 2010) is a global outlook on web mapping solutions (both client and server).
comparison of FOSS libraries by German Carrillo can be found here:

The list of libraries so far:

Google Maps
Microsoft Virtual Earth
MapQuest
Leaflet - "The comment smaller, faster, newer, and more straightforward can also be read as less features and less tested." -Geographika (see below)
ArcGIS API for JavaScript - Works best with ArcGIS Server (see below). Google Maps and Bing maps extensions are also available, letting you use the ESRI API with Google/Bing maps (though this is true of most libraries).
Yahoo Map API's
Via Michelin
OpenLayers - Extensive documentation and a good amount of functionality plus the ability to use different map providers.
Mapquery - MapQuery has been released and now has some useful documentation It has the very worthwhile goal of combining OpenLayers and jQuery. If you're particularly keen on the idea of OpenLayers + jQuery, or if you want to contribute to a JavaScript Mapping Library then get involved and contribute your efforts. However if you just want to be an end user, or are new to this area it may not be for you.
Mapstraction - Makes things very simple, especially working with multiple basemap providers. However it is still a work in progress and the functionality is lacking in places, as is the documentation. (E.g. "A GeoJSON object with the type “FeatureCollection” is a feature collection object." Not very informative.) It appears that it is still being actively developed but as of 4/4/11 there hasn't been a commit on Github since January.
deCarta - Has a mobile and desktop javascript - first is HTML5/CSS3 compliant and the second has more browser compatibility. Source code provided. Friendliest developer terms for a commercial API. You are allowed to brand the map and there are several different map styles. You can choose NAVTEQ or OSM data. They also have several Mobile APIs as well. - edited by TheSteve0 - a deCarta employee
Cloudmade
Polymaps - Makes it very easy to composite raster and vector data from many different sources. Lets you easily add your own colouring, grouping, and interaction. Runs quickly, manages background tile loading well, and is only 30k of Javascript. One potential down side: it uses SVG which means it does not and will not work in MSIE 7 or 8. It works great in every other browser and should work in IE9
Jump - jump is a light weight maps library that works on its own, meaning, it is not a wrapper for OpenLayers or GoogleMaps API. Currently it is under development, but a lot of essential features work well.
ModestMaps - Another smaller, faster, newer JS mapping library from the makers of Mapbox and TileMill.
Mapiator

OpenLayers is the one I am currently using. You can do a lot with it and it supports most data types. However it's not the best for everything. For example, Leaflet seems smoother in many ways, with image fading and other visual tweaks. If you're into jQuery you might like to check out MapQuery which is like a combination of jQuery and OpenLayers.


texts:
I want to get the native SRS/EPSG-Code of a layer deployed in GeoServer. Maybe similar to this solution: https://stackoverflow.com/questions/32766126/get-the-epsg-code-from-openlayers-3-layer but with Openlayers 2 or the RESTful API of GeoServer.

Is there a way to achieve this?
GeoServer itself displays those native EPSG codes when you look at your layers on the right hand side.

So somehow this must be possible but I just can't find the answer.


texts: Load multiple files from folder for further processing

texts:
I'm working on a processing model and I'm looking for a possibility to load multiple shape-files from an Input-folder for further processing (first clip them all to a certain extend, then assign them to variables for further processes etc.) with a python script.

The problem is that the number of files in the Input-folder can vary from model-run to model-run. So far, the user has to define a certain number of shapefiles as Input.

Any Ideas about that?


texts: How to enclose clustered points with polygon?

texts:
I have a raster with forest (1), no forest (2) data. "1" is represent single/few trees and there is gaps between them. I want to create a polygon around these trees to a forest map. It is not a dense forest, so it is not possible to filter the data to fill in the gap between trees.
Is there any tools in ArcMap 10.0 or extensions that could work to use? I have seen a tool doing this a while ago, don't remember if it was a ArcMap tool/extension or if it was in another GIS software.

Regards


texts: Geomorphons for a geomoron: running grass command in external python

texts:
I'm a fairly competent Linux and Python user and looking to utilize the GrassGIS geomorphon function (https://grass.osgeo.org/grass78/manuals/r.geomorphon.html) to parse out subaqueous features.  I've never used GRASS and don't quite understand the architecture and work flow, so please bear with me.

In an ideal world I install the GRASS backend only via pip eg pip install pyGrass and call import pyGrass and  pyGrass.r.geomorphon(arg1, arg2) (or something similar) similar to an interaction like is seen with the C based opencv library (https://opencv.org/) eg. pip install opencv then import cv2.  My understanding is that not possible with the GRASS-GIS as it stands.

For example my envisioned workflow:

import pickle
from matplotlib import pyplot as plt
url= "https://chldata.erdc.dren.mil/thredds/dodsC/frf/geomorphology/DEMs/surveyDEM/FRF_20191114_1177_FRF_NAVD88_LARC_GPS_UTC_v20191115_grid_latlon.nc"
ncfile = nc.Dataset(url)
lon = ncfile['longitude'][:]
lat = ncfile['latitude'][:]
elev = ncfile['elevation'][:].squeeze()

returnedArguments = pyGrass.r.geomorphon(elevation, arg1, arg2)
plt.figure()
plt.pcolormesh(lon, lat, elevation)
for idxType in returnedArguments:
    plt.pcolormesh(lon, lat, featureType[idxType])

I've found a few resources that I've begun cobbling together to work towards my end, none of which seem to illustrate my goals directly.  It looks like most everything in the documentation is Python2 and I'm using Python3, CentOS

import os, sys, subprocess, pickle
url= "https://chldata.erdc.dren.mil/thredds/dodsC/frf/geomorphology/DEMs/surveyDEM/FRF_20191206_1179_FRF_NAVD88_LARC_GPS_UTC_v20191209_grid_latlon.nc"
ncfile = nc.Dataset(url)
lon = ncfile['longitude'][:]
lat = ncfile['latitude'][:]
elev = ncfile['elevation'][:].squeeze()

location_path = os.path.join(gisdb, location)
if os.path.exists(location_path):
    shutil.rmtree(location_path)
# Create new location from EPSG code, and exit
startcmd = ("%s -c %s -e %s" % (grass7bin, epsg_code, location_path))
print(startcmd)
p = subprocess.Popen(startcmd, shell=True,
                     stdout=subprocess.PIPE, stderr=subprocess.PIPE)
###########################
grassPath = "/home/spike/grassStuff"
#################################################
# from https://github.com/wenzeslaus/try-grass-in-jupyter/blob/master/notebook.ipynb
# create GRASS GIS session
gisbase = str(subprocess.check_output(["grass78", "--config", "path"]), 'utf-8').strip().split()[-1]
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))
# import GRASS GIS packages we need
import grass.script as gs
import grass.script.setup as gsetup
import grass.script.array as garray
rcfile = gsetup.init(gisbase, grassPath, location='demo')
# we want functions to raise exceptions and see standard output of the modules
gs.set_raise_on_error(True)
gs.set_capture_stderr(True)
# simply overwrite existing maps like we overwrite Python variable values
os.environ['GRASS_OVERWRITE'] = '1'
# from https://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library#Interfacing_with_NumPy
a = garray.array()   # I error out here, but this is not the focus of my question, i'm curious why
a.read(elevation)

As shown above, I would love to get here and is the main focus of my question:

# I'm guessing that I have to write the data to a raster
gs.grass.run_command(gscore.g.region raster=elevation -p)
# then gs.execute or some equivalent that returns data
returnedArgs = grass.run_command(r.geomorphon elevation=elevation forms='all')
# then return values back as numpy array
dataIAmInterestedIn = a.write(returnedArgs)

but I'm quite sure that's totally wrong. What should I do instead?


texts:
I've never used TransCAD before. I have what I believe to be a TransCAD matrix file called data.mtx which I would like to be able to convert to some other file formats. I would prefer to do this at the Unix command line. I could use TransCAD directly to convert this to another format, but I'd like to be able to do this in an automated way via a shell script. I tried to search for TransCAD file formats and came across this page:

Map Your Own Data with TransCAD

Unfortunately I didn't see any occurrences of the strings 'matrix' or 'mtx'. Here is my first attempt at trying to understand what this file is:

$ file AM.mtx
AM.mtx: ARC archive data, squeezed

$ nomarch AM.mtx
nomarch: error reading data (hit EOF)

Some other stackexchange posts mentioned that the FreeARC utility may succeed in extracting an ARC file where nomarch fails, but their download links seem to be dead.

So what's the deal? Is this not an ARC file? Is there some documentation on the TransCAD matrix format? Are there any open-source command-line tools for working with TransCAD files?


texts: Shortcuts:run layer actions without click

texts:
I have added an action on a layer(Layer_1) and assigned a shortcut on the keyboard to execute this action using 'Keyboard shortcut settings ..' from the toolbar, however, by pressing the number '2' which is the shortcut used by the user keyboard only activates the actions icon but the script associated with that action is not executed. To run it is necessary to click on any feature of the active layer.This is the script that must be executed:

from qgis.utils import iface
ly2 = QgsMapLayerRegistry.instance().mapLayersByName('Layer_2')[0]
iface.setActiveLayer(ly2)

Would it be possible to run this script without having to click on any feature of the active layer with that action?


texts: Scaling attributes when clipping / intersecting in QGIS?

texts:
I'm working on a project involving a ratio of frequent travellers / population divided in Forward Sortation Areas from Canada Post (FSAs are the first three characters of a postal code: A1A).

I created an 8km buffer around stations and then Intersected the FSA data with the 8km buffer, only to realize that the attributes were completely imported.

It seems to me that there is an option in ArcMap that scales the attributes proportionaly to their size when clipped, is there such a thing in QGIS or do I have to do this manually?


texts: sld density heatmap with geoserver without weight atribute

texts:
I try to generate a density heatmap to display my points density but I don't know what to write in the weightAttr area of the tamplate furnished in the doc because I just have plenty of points but they don't have weight attribute.

I can generate it with QGIS but the sld export isn't available for heatmap.

I can also generate it using openlayers from geojson file but I have more than 1.000.000 points and it crashes.

I think it's better to use postgis and geoserver to generate wms to display so many points.


texts:
I am new to using spatial data in R. So far the tutorials I have seen consist of using maptools, map, ggplot2 and rgdal to open ".shp" files along with ".csv" files to plot the polygons/ points using latitude and longitude.

I have always used ArcGIS and presently have these set of files: US-101.shp, US-101.shx, US-101.dbf and US-101.mxd

I would like to open these files so that I could use the shapefile along with a data frame related to it. I can open the shapefile and also see the table which has "ID", "LAYER", "COLOR", "LINETYPE", "LINEWIDTH" and "ELEVATION" variables.

I am not sure how to use it to plot the shapefile. In ArcGIS the layer is easily visible once I open .mxd file.


texts:
I need help changing this code to do:

1) Within a specified directory, reproject all .tif and/or .shp files to:  NAD_1983_UTM_Zone_16N without needing to input what the old projection is.

2) If already projected correctly, skip.

I attempted to search and manipulate existing code already out there, and found this:
How do I use python to batch process the reprojection of a group of shapefiles?

import sys, arcpy, os

mypath = r"C:\Users\Robby\Desktop\Test"
shapefiles = []
arcpy.env.workspace = mypath

for ThisFile in os.listdir(mypath):
    thisName,thisExt = os.path.splitext(ThisFile)
    if thisExt.upper() == ".TIF": # or thisExt.upper() == ".SHP": # include your shape files as well
        shapefiles.append(ThisFile)

out_workspace = r"C:\Users\Robby\Desktop\Test"
out_cs = ''
template = "C:\\"
transformation = "NAD_1983_UTM_Zone_16N" # that's not a transform, but I'm not an  expert on North American projections
OutSR = arcpy.SpatialReference(26916)    # http://spatialreference.org/ref/epsg/nad83-utm-zone-16n/, EPSG:26916

for ThisFile in shapefiles:
    desc = arcpy.Describe(ThisFile)
    SR = desc.spatialReference
    if SR.factoryCode != 26916: # http://spatialreference.org/ref/epsg/nad83-utm-zone-16n/, EPSG:26916
        arcpy.Project_management(ThisFile,out_workspace + "\\" ThisFile,OutSR)


texts: Shapely command error: from shapely.geometry import LineString

texts:
I am trying to run a program but getting the following error when I am using the command:
from shapely.geometry import LineString

Error: FileNotFoundError: Could not find module
'C:\Users\SWWB\Anaconda\Library\bin\geos_c.dll' (or one of its
dependencies). Try using the full path with constructor syntax.

Can anyone help me with this?


texts: Conditional spatial join based on largest overlapped area

texts:
I have two layers (Buildings and Census areas) which are overlapped. I want to join them with a spatial join but, in some cases, buildings overlap two or more census units and I don't know which census unit associating to. So I thought that I could associate the building to the census area that is overlapped for the largest part. For doing this, I first want to evaluate if the intersected area is higher than the half of the building area. I tried in two ways (see code below) but they return different results.

does someone know why the two codes return different results?

  1st trial

featuresC= Census_layer.getFeatures()
featuresB=Building_layer.getFeatures()
count=0
for fbuild in featuresB:
    BuildingGeometry=fbuild.geometry()
    if intersection.area()>0.5*BuildingGeometry.area():
        print intersection.area(), "featInters"
    else:
        print BuildingGeometry.area(), "featBuild"

  2nd trial

featuresC= Census_layer.getFeatures()
featuresB=Building_layer.getFeatures()
count=0
for fbuild in featuresB:
    BuildingGeometry=fbuild.geometry()
    if intersection.area()*2.0>BuildingGeometry.area():
        print intersection.area(), "featInters"
    else:
        print BuildingGeometry.area(), "featBuild"

and, for instance, the results are:

         **1st trial             2nd trial**

    717.0162316 featBuild   717.0162316 featBuild
    215.4575371 featBuild   215.4575371 featBuild
    78.8431039  featBuild   50.23154609 featInters


texts: add geometry to MULTILINESTRING in cartodb

texts:
I have a table that contains a column with MULTILINESTRINGS (which are the coordinates from two places (i.e. long0,lat0,long1,lat1), like this:

MULTILINESTRING((97.0253024 4.3685491,8.5618685 47.4524688))

When I now click on "Map View", nothing gets displayed because "the_geom" in my table is null. My question now is: how can I add geometry information to the column with the MULTILINESTRINGs so that they get displayed on the map? I tried by updating "the_geom" like this:

UPDATE table_name
    SET the_geom =
      ST_SetSRID(column_with_multilines,4326)

But I only get "undefined rows affected". Could anyone please help me?


texts: Does SQL Server need extensions to do spatial?

texts:
Does SQL Server need extensions to do spatial like Oracle does?
Just wondering if there will be extra costs for implementing a SQL Server backend to SDE.


texts: Is there a way to create contour map from a given shaded relief map?


texts:
I am trying to mask out a section of an image I created with an algorithm. I create a geometry of the section I want to mask out and clipped that with an image of 1's to create my mask image. But when I attempt to mask my original image it only clips the part I'm trying to remove rather than masking out the section and retaining the rest.

how can I remove the pixels from my mask and retain the rest of the image?

link: https://code.earthengine.google.com/844c40e72f98761c0bbe46a003989213


texts: Geometry not showing when importing Dutch BRT-files (GML)

texts:
After importing a couple of gml files (Dutch BRT-files) into QGIS I noticed that some of the layers didn't seem to contain any geometry, the database (attribute table) is there, but no geometry is showing on my screen.

When I import the gml separately, it shows up just fine.It looks like the gml files are interfering with each other?

Any thoughts?


texts: Extract latitude and longitude from rotated coordinate system

texts:
I'm new to working with geographic data so please bear with me.

I have the following data:

bnds,projection_x_coordinate,projection_y_coordinate,realization,fog_area_fraction,lambert_azimuthal_equal_area,projection_y_coordinate_bnds,projection_x_coordinate_bnds,forecast_period,forecast_reference_time,time
0,-1158000.0,-1036000.0,12,0.0,-2147483647,-1037000.0,-1159000.0,1 days 06:00:00.000000000,2019-03-22 19:00:00,2019-03-24 01:00:00
0,-1158000.0,-1036000.0,13,0.0,-2147483647,-1037000.0,-1159000.0,1 days 06:00:00.000000000,2019-03-22 19:00:00,2019-03-24 01:00:00
0,-1158000.0,-1036000.0,14,0.0,-2147483647,-1037000.0,-1159000.0,1 days 06:00:00.000000000,2019-03-22 19:00:00,2019-03-24 01:00:00
0,-1158000.0,-1034000.0,12,0.0,-2147483647,-1035000.0,-1159000.0,1 days 06:00:00.000000000,2019-03-22 19:00:00,2019-03-24 01:00:00

from which I would like to determine the appropriate lat/long coordinates. Note that the data is in the RotatedGeogCS format and I have x- and y-projections along with a data value for "Lambert Azimuthal Equal Area" projection.

Does anyone know how to convert from the given data to what I want? Is that even the appropriate question to ask here? The end goal is to present the fog_area_fraction data within the region it has been calculated for.

For those curious this data was extracted from the MOGREPS-UK dataset from AWS, processed using libraries for netCDF, xarray and Iris and then dumped to CSV.

EDIT: Here's a sample of the cube in the original form:

    Dimensions:                       (bnds: 2, projection_x_coordinate: 1042, projection_y_coordinate: 970, realization: 3)
Coordinates:
  * realization                   (realization) int32 12 13 14
  * projection_y_coordinate       (projection_y_coordinate) float32 -1036000.0 ... 902000.0
  * projection_x_coordinate       (projection_x_coordinate) float32 -1158000.0 ... 924000.0
    forecast_period               timedelta64[ns] ...
    forecast_reference_time       datetime64[ns] ...
    time                          datetime64[ns] ...
Dimensions without coordinates: bnds
Data variables:
    fog_area_fraction             (realization, projection_y_coordinate, projection_x_coordinate) float32 ...
    lambert_azimuthal_equal_area  int32 ...
    projection_y_coordinate_bnds  (projection_y_coordinate, bnds) float32 ...
    projection_x_coordinate_bnds  (projection_x_coordinate, bnds) float32 ...
Attributes:
    history:                      2019-03-22T22:54:56Z: StaGE Decoupler
    institution:                  Met Office
    mosg__forecast_run_duration:  PT126H
    mosg__grid_domain:            uk_extended
    mosg__grid_type:              standard
    mosg__grid_version:           2.1.0
    mosg__model_configuration:    uk_ens
    source:                       Met Office Unified Model
    title:                        MOGREPS-UK Model Forecast on UK 2 km Standa...
    um_version:                   11.1
    Conventions:                  CF-1.5, UKMO-1.0


texts:
I have several dozen points plotted across the UK. Each point represents a facility that serves local people. To estimate the catchment areas of each facility I created Voronoi polygons around them, and clipped the polygons to the UK borders (using a UK shapefile). I then coloured each polygon according to an attribute of the point data.

A cropped part of the map is here: (Wales and Devon/Cornwall)

Is there any way to make the Voronoi respect the natural borders of the UK (i.e. limit them to my shapefile during Voronoi creation)? The white on my map is the sea, and some polygons have crossed the Bristol Channel. This is implausible for a catchment area.

Is this possible, or is there any alternative? I have used QGIS so far, but am open to any option, code or GUI, preferably free.

Edit: There are some questions Constructing a Voronoi diagram using a complicated travel time metric and Is there a described GIS method for creating catchment areas based on nearest driving distance? that talk about building catchment areas from driving times. That is a possibility that I am considering if this is not possible, but a Voronoi algorithm that respects borders would be simpler.


texts:
I'm trying discover how to read the time:unit variable information information from a netCDF file and have run out of ideas.  Without this info i cannot auto-magically process the files by date (which is my objective).

Any suggestions?

I am using the latest and greatest GDAL with C# bindings (gdal_csharp, runtime version v4.0.30319, version 1.0.5507.11621).  VS2010 64 bit.


texts:
I'm trying to zoom to a feature just like in Arcmap but in VB.net.
I'm using Axmapcontrol to display the map.
I've tried the following:

Dim pCmdItem As ICommandItem
Dim pUID As New UID
pUID.Value = "{AB073B49-DE5E-11D1-AA80-00C04FA37860}"
pCmdItem = CommandBars.Find(pUID)
pCmdItem.Execute()

but it's showing the following error:

Reference to a non-shared member requires an object reference.

I've found a solution that suggets to change mytype from empty to WindowsForms in the vbproj file, but it already contains the correct value.
Is there any way to zoom to a feature in an Axmapcontrol in VB.net?


texts:
ERROR 3: Free disk space available is 55518121984 bytes, whereas 69156878736 are at least necessary. You can disable this check by defining the CHECK_DISK_FREE_SPACE configuration option to FALSE.

I receive this error when trying to re project a .tif file on my EC2 instance. I understand that warping a file requires more memory than the size of the file itself. Is there any way to decrease the amount of memory required?

Here is my command:

gdalwarp --config GDAL_CACHEMAX 2048 -wm 512 -wo NUM_THREADS=ALL_CPUS -srcnodata nodata_masking_values -dstnodata 0 input_file.tif output_file.tif


texts: How to tie an ArcGIS Server QueryTask result back to the calling layer?

texts:
Using ArcGIS Server JS API 4.x I'm running a series of QueryTasks in succession, and need to tie the result of the query back to the calling layer. However, since the QueryTasks execute asynchronously I'm not sure which result corresponds to which QueryTask.
To illustrate the issue I created a simple example on CodePen which writes the results of the query to the console.
Open the debugging tools and note that the result of the QueryTask has no information as to which layer the results apply to:
require([
  "esri/tasks/QueryTask",
  "esri/tasks/support/Query"
], function(QueryTask, Query) {

 var config = {
   "layer1": {
        "url": "https://sampleserver6.arcgisonline.com/arcgis/rest/services/WorldTimeZones/MapServer/0"
    },
    "layer2": {
        "url": "https://sampleserver6.arcgisonline.com/arcgis/rest/services/911CallsHotspot/MapServer/1"
    },
   "layer3": {
        "url": "https://sampleserver6.arcgisonline.com/arcgis/rest/services/LocalGovernment/SewerManholeInspections/FeatureServer/0"
   }
}

for (layerDef in config){
    var queryTask = new QueryTask({url: config[layerDef].url});
    var query = new Query();
    query.where = "1=1";

    // These execute asynchronously so there is no way to tie the result back to the calling query
    console.log("about to check count for", layerDef);

    queryTask.execute(query).then(function(results){
        // How to detect which layer this result corresponds to?
        console.log("Count:", results)
    });
}

How can I determine to which layer the results correspond? I guess I need a separate callback function for each queryTask - but I can't see how to specify this.
PS. Here I'm using queryTask.executeForCount but the same applies if I use a simple queryTask.execute, and the principle should hopefully be the same anyway.
PPS. An option could be to make an AJAX request directly against the server URL, but due to the environment I'm working in that fails due to CORS issues, so I'm stuck using the QueryTask approach


texts: Download Polygon/Administrative Boundaries using WikiMapia

texts:
Is there any way to download administrative boundaries / polygons for a small city through Wikimapia ??
i tried to download the boundaries but the results were points not polygons.
i need to download the polygons from wikimapia


texts: Is there an Attribute Table Feature limit?

texts: I am copying a table with 11,000+ features but only about 700 show up on the QGIS attribute table. Is there a limit? how can I avoid this problem?


texts: Solving error 0x80131523: Entry point was not found activating SOE

texts:
I have ArcGIS Server Object Extension (SOE) written in .NET, C#. Target framework is .NET 4.5.

Extension works OK at many customers but there is an issue at one customer who has class extensions targeting .NET 2.0.
Uploading SOE file to ArcGIS Server - OK.
Enable our extension for map service and starting it, an error occurs:

  Failed to initialize server object 'MapServiceName': 0x80131523: Entry
  point was not found.

Not much information what went wrong... I tried to turn log level to debug - ArcGIS doesn't show more details.

My question is how to get more information about such error?

Turn on catalina error log for ArcGIS internal web server? Other methods?

The issue is not the mix of .NET 2 and .NET 4 frameworks. The previous SOE version was also targeting .NET 4.5 but it worked side by side with these .NET 2.0 extensions fine.

HRESULT 0x80131523 is COR_E_ENTRYPOINTNOTFOUND - Could not find the specified DllImport entrypoint.

Created REST SOE proxy to invoke AssemblyLoad.Load() inside ArcSOC and tried to load real SOE dlls. All were loaded OK.

Removed 3rd party extensions while found the single extension showing this behavior. It is ArcFM Solution Object Reader 10.2.1c for ArcGIS Server (by Schneider Electric).


texts: Renaming shapefiles using Windows Explorer?

texts:
I need to rename a shapefile.

Is it ok to simply use Windows Explorer to manually rename each of the three related files that comprise the shapefile?

I am using QGIS.


texts:
Is it possible to get the name of an area providing just the bounding box?

I'm not interested in specific points, what I'm trying to figure out is to get the name of the maximum polygon (I'm thinking in OSM here) within the bounding box I already have.

Clarifying example: Send bounding box containing part of New York and get "Brooklyn" as a response.


texts: How do you Aggregate a set of points that lie within different sections of a polygon?

texts:
I require some help with an ArcGIS 10.2.1 map that I am trying to create. I want to aggregate a set of points into bigger points depending on the number of those points that are inside a section of a polygon.

I have a shapefile with a set of points that represent service requests in New York City. I also have a shapefile that is a polygon of New York City divided up into different census sections. I would like to combine the individual points into larger points in relation to how many of those points are within a specified section of the polygon. I want to do this for each individual section of the polygon.

So that the resulting map will have one point for each section of the polygon. The size of the point will differ based on how many points were aggregated together to create that resulting point.


texts:
I would like to read multiple features with ol.format.WKT.
In the documentation, it says that the function readFeatures takes an array as parameter. I tried the following:

let wkt = new ol.format.WKT;
wkt.readFeatures(["POINT(-378147.31208516686 599021.3516006605)", "POINT(-378300.01500700554 599159.7386235767)"]);

But I get an error :

Uncaught Error: Unexpected `` at position 0 in ``
    at Tu (ol.js:659)
    at Ku.k.Gd (ol.js:655)
    at Ku.k.Ng (ol.js:655)
    at Ku.k.Qa (ol.js:562)
    at <anonymous>:1:5

Has anyone any clue ? I tried searching on the internet but no exemples are available for this specific method, and there are none either in the documentation...

EDIT
The function doesn't takes an Array as parameter but a string or an Object (if I red well). What are the syntaxes to use those inputs ? I had no success so far with GEOMETRYCOLLECTION string (see below).


texts: Can multiple point symbols be stored in one QGIS layer?

texts:
I'm new to QGIS, indeed fairly new to GIS systems in general.
I have been asked to create about 20 small maps for a local historical society.  The maps are for walkers to follow various routes passing locations of historic interest.

Each of the maps is to show numbered "points of interest" (POIs) - some maps have only 5, others have 12 or more.  I have found and imported the symbols I want into QGIS.  Each of the symbols consists of an integer enclosed in a circle.

My first attempt at getting the symbols onto the maps was to create a separate layer for each symbol.  This seems to work OK but means that a large number of layers have to be created for each map - one for each POI.
Is there a more efficient way of doing the job?

The POIs have only the 2 attributes, ID and the relevant integer.  I was hoping to be able to just drag them to their correct location on the map.

I have SVG images of numbers 1 to 15 each in a circle.  I want to be able to place them on my map at specific locations to indicate points of interest.  How do I do this please?


texts: Why is this happening? I show GeoJSON some coordinates with OpenLayers by request Ajax

texts:
I want to display a GeoJSON on OpenLayers with a request from Ajax from the database:

I store in the geojson variable all the information and I icercated in two ways:
1. Directly from the already read variable -> in this case it works
2. From the external source in the json / geojson.json folder -> that does not work

<script type="text/javascript">

      proj4.defs('EPSG:31700','+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 +y_0=500000 +ellps=krass +towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs');
      ol.proj.proj4.register(proj4);

var geojson = $.getJSON('json/geojson.json', function(data) { });
var geojson = {"type":"FeatureCollection","features":[{"type":"Feature","id":"2","geometry":{"type":"LineString","coordinates":[[698282.236,377375.758],[698289.591,377372.136],[698183.88,377182.298],[698182.72,377197.045],[698282.236,377375.758]]},"crs":{"type":"EPSG","properties":{"code":"31700"}},"geometry_name":"wkb_geometry","properties":{"ogc_fid":"2","area_id":"2","building_id":"2","type":"test","area":"test","owner":"test","description":"test","administration":"test","address":"test","property":"test","layer":"Green Area","extendedentity":null,"linetype":"test","entityhandle":null,"text":"test","wkb_geometry":"01020000000500000000FED478544F2541002731087F081741408D972E634F25418195438B70081741295C8FC28F4E254181E92631790517410AD7A3708D4E2541027B142EB405174100FED478544F2541002731087F081741"}}]}

console.log(geojson);

var features = (new ol.format.GeoJSON()).readFeatures(geojson);

features.forEach(function(feature){
    feature.setId(undefined);
    feature.getGeometry().transform('EPSG:31700', 'EPSG:3857');
});

var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM(),
        }),
        new ol.layer.Vector({
            source: new ol.source.Vector({
                features: features
            })
        })
    ],
    target: 'map',
    controls: ol.control.defaults({
            attributionOptions:
                ({
                    collapsible: false
                })
        }).extend([
            new ol.control.OverviewMap(),
            new ol.control.FullScreen(),
            new ol.control.ScaleLine()

        ]),
        interactions: [
            new ol.interaction.DragRotateAndZoom(),
            new ol.interaction.MouseWheelZoom(),
            new ol.interaction.DragPan(),
            new ol.interaction.DoubleClickZoom()
        ],
    view: new ol.View({
        center: ol.proj.transform([698282.236, 377375.758], 'EPSG:31700', 'EPSG:3857'),
        zoom: 14,
    })
});

    </script>

The GeoJSON content is the same does not differ.


texts: Difference between latitude and longitude

texts:
It occurred to me that latitude and longitude are not homologous in that the lines of latitude lie on parallel planes or slices through Earth, while lines of longitude are not parallel and meet at the poles, forming segments through Earth.  Obviously the Earth's rotation gives us two points that make sense to use as intersection points, but it struck me that there could be potential advantages for using arbitrary points for intersection points for longitude, or for disregarding the poles for latitude.   Not knowing the names of these I was not able to find any information online.

Is there a name for this difference, or names for the two types?  Are there any coordinate systems that use parallel longitudes (or non-parallel latitudes) and what would the implications of this be?


texts: Calculating green metropolitan urban areas

texts:
I need to calculate green spaces in the metropolitan area for 155 cities around the world in QGIS.

Greenery would include everything from parks, trees, playgrounds, golf, vineyards, etc. Any advice where to start?

I have found a shapefile from NYU for metropolitan areas and now looking for each city under their departments, but already see an issue with that.

I wanted to do data scraping from OSM maps (it would be not all including). Also tried iTree even for approximate results but there were cities with missing or unclear satellite images.


texts: Passing location data to a global variable in JavaScript

texts:
I am trying to pass the lat and lon of a user location to a global variable within JavaScript.
The variable is instantiated at the top of the script and given value within the location function.
//lat long variables of user
var lon, lat;

//User location
var Position,
  Accuracy,
  test = null;
function onLocationFound(e) {
  if (Position) {
    pub_map.removeLayer(Position);
    pub_map.removeLayer(Accuracy);
  }
  test = e.latlng;
  console.log("user lat long in script is: " + test);
  lon = test.lng;
  lat = test.lat;
  var radius = e.accuracy / 2;
  //pop up stating the users position
  Position = L.marker(e.latlng)
    .addTo(pub_map)
    .bindPopup(
      "
You are here

You are  within " +
        radius +
        " metres from here
"
    )
    .openPopup();
  Accuracy = L.circle(e.latlng, radius).addTo(pub_map);
  return lon;
  return lat;
}

//error message if user location cannot be found
function onLocationError(e) {
  alert(e.message);
}
//display the layers to the map
pub_map.on("locationfound", onLocationFound);
pub_map.on("locationerror", onLocationError);

//locate function
pub_map.locate({ setView: true, maxZoom: 16 });

console.log("user lat is: " + lat);
console.log("user lon is: " + lon);

However when trying to call the value outside the script the value is undefined. I know the variable is global and should be accessible by all aspects of the script but for some reason when logging the data outside it still says undefined.
I have tried wrapping the whole script within a function (a bodge if I do say so) and this did not work.
I have tried calling it as a reference but this too did not work.
Any ideas or reasons to why this is and how to fix it?
The image below show the logs as you can see inside the function the value is obtained but it is not passed outside of the function.


texts: Assign RasterB values to RasterA where B != no value in QGIS Raster Calculator

texts:
QGIS raster calculator is confusing me and I have no idea how to perform a very simple task.

I have 2 rasters, A and B:

Raster A has 9999 values on all pixels.
Raster B has different values including no data pixels.

I want Raster A pixel to = Raster B pixel where B != null

Which means Raster C(output) is 9999 Except where Raster B has values


texts:
Have some coordinates for a spdf of the Green Mountain National Forest in Vermont in R. Downloaded from Vermont.gov website. http://geodata.vermont.gov/datasets/vt-green-mountain-national-forest-mgmt-areas-2006-polygons?geometry=-74.248%2C43.079%2C-71.347%2C43.777
The coordinates look this:

 gmnf <- matrix(c(437238.82, 26905.22, 514258.28, 184930.56), nrow = 2, dimnames = list(c("x","y"),c("min","max")))

I thought they were in DMS and I put them in a few online calculators but the output values don't look right. They should fall within the range of Vermont but they don't. From a spdf for Vermont:

vt <- matrix(c(-73.43774, 42.72685, -71.46455, 45.01666), nrow = 2, dimnames = list(c("x","y"),c("min","max")))

Based on the metadata (linked to on the above website) these are coordinates in decimal degrees. When I change them and try to plot(gmnf)nothing happens. Tried with and without changing the crs but still nothing.

gmnf@bbox <- matrix(c(-73.454162, 42.722789, -71.465281, 45.018361), nrow = 2, dimnames = list(c("x","y"),c("min","max")))
gmnf@proj4string <- vt@proj4string
plot(gmnf)

What types of coordinates is gmnf in and how do I change it to decimal degrees?

Dropbox folder with the spdf for gmnf and vt:
https://www.dropbox.com/sh/i94ondu3cihglef/AAA6jZV2cwpersPNs_9oUP6ea?dl=0

Update: Solution is painfully simple. Noob

gmnf <- readOGR(dsn = "gmnf_shape", layer = "VT_Green_Mountain_National_Forest_Mgmt_Areas_2006__polygons")
latlong = "+init=epsg:4326"
gmnf <- spTransform(gmnf, latlong)


texts: Why isn't Mapbox viz styled how it is in projects?

texts:
When using the Mapbox "Projects" editor, I'm finding the previews do not style the same as the end results on the page. Is anyone experiencing similarly?

In Projects, my map looks like on the left. On the right is how the project appears when visualized:

Interestingly, on the Projects page, the preview looks like the visualization on the right. It appears like this:

Does anyone know what's happening?


texts:
I have population distribution data as a raster layer and exposure contours as vector (polygon) layers. What I'd like to do is use the vector layer to estimate the population indicated from the raster layer that fall underneath the polygon.

Can someone shed some light on this?

If it isn't straight forward, I am not too bad with Python so I'd be more than comfortable writing something that does this.


texts: 80% of data has a spatial component, says who?

texts: This quote about how 80% of all data has a spatial component has been bandied around everywhere. Introduction to GIS classes, product descriptions, talks, etc. Would anybody know who first said it? Or  A link to original article would be nice.


texts: How do I generate 3d WKT FROM A Sql Server geometry or geography types

texts:
If I have the point POINT (-75.45454 25.12345 1024 ) in a geometry or geography column in SQL Server, but render it with STAsText() I lose the Z dimension.

SELECT geography::STGeomFromText('POINT (-75.45454 25.12345 1024 )', 4326).STAsText() AS point
POINT (-75.45454 25.12345)

How do I get the full WKT? There is nothing listed in OGC Methods on Geography Instances.


texts: Numbering of objects based on several fields in ArcGIS with ArcPy

texts:
I have a feature layer. I need to number objects within the numbering of houses on each street. The first field is the city code, the second field is the street code, the third - house number. A little confused how can I number this. The picture shows how the final "VOPROS2" field should appear.

import arcpy
numarr=[]
curs=arcpy.da.SearchCursor("Lot",["SOATO","Address3","adr2","OBJECTID"])
for row in curs:
    numarr.append(row)
numarr.sort(key=lambda x: (x[0],x[1],x[2],x[3]))


texts: ST_Intersects with subquery in Spatialite

texts:
I am having issues with what seems to be a very simple spatial query with a sub-query in Spatialite. I am trying to find all the buildings (poly) that intersect with a service area (poly), but only buildings that have a specific subtype code and with stories greater than a value (see query below). I am assuming it is a syntax error, but cannot seem to figure it.

SELECT b.pkid, b.subtype, b.stories
FROM bldg AS b, service_area AS a
  (SELECT b.subtype, b.stories
   FROM bldg
   WHERE subtype IN ('x', 'y', 'z') AND stories > 16)
WHERE ST_Intersects(b.geometry, a.geometry)

ERROR

Query execution error on CREATE TEMP VIEW _tview AS SELECT b.pkid, b.subtype, b.stories
FROM bldg AS b, service_area AS a
  (SELECT b.subtype, b.stories
   FROM bldg
   WHERE subtype IN ('x', 'y', 'z') AND stories > 16)
WHERE ST_Intersects(b.geometry, a.geometry): 1 - near "(": syntax error

EDIT

Query to select desired subtype and buildings > 16 stories

ST_Intersects query - buildings and service area (buildings shorter than 16)

I did get the following query to work, but its not a subquery and not the most optimal:

 SELECT b.pkid, b.subtype, b.stories
 FROM bldg AS b, service_area AS a
 WHERE ST_Intersects(b.geometry, a.geometry) AND b.subtype IN ('x', 'y', 'z') AND b.stories > 16


texts: Wedge of missing data reprojecting to stereographic

texts:
I am trying to reproject an equirectangular projection image into stereographic around the north pole of mars using gdalwarp and I am not sure what is happening in my output image.

The input image is a normal looking rectangular image with no black regions.  Zooming in to the missing wedge region I can see that not all of the pixels from the input image are present in the output image.

My gdalwarp command:

gdalwarp map_upper_crop.tif map_north_pole.tif -t_srs "+proj=stere +lon_0=0 +lat_0=90 +a=3396200 +b=3376200 units=m"

gdalinfo output on the input image:

Driver: GTiff/GeoTIFF
Files: projection_space_basemap_top.tif
Size is 11520, 1035
Coordinate System is:

PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",3396200,169.8100000000004]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

PROJ.4 string is:
'+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396200 +b=3376200 +units=m +no_defs '
Origin = (-10669477.099999999627471,5334738.599999999627471)
Pixel Size = (1852.339774305555466,-1852.339791666666542)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=1024
  TIFFTAG_YRESOLUTION=1024
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-10669477.100, 5334738.600) (179d59'59.99"E, 90d 0' 0.01"N)
Lower Left  (-10669477.100, 3417566.916) (179d59'59.99"E, 57d39'22.50"N)
Upper Right (10669477.100, 5334738.600) (179d59'59.99"W, 90d 0' 0.01"N)
Lower Right (10669477.100, 3417566.916) (179d59'59.99"W, 57d39'22.50"N)
Center      (       0.000, 4376152.758) (  0d 0' 0.01"E, 73d49'41.26"N)

Any idea what is going wrong here?


texts: QGIS - Topological colour doesn't add new layer

texts:
I have been reading the "Handling Invalid Geometries" topic in the QGIS Tutorials & Tips section, to which the link you can find here: https://www.qgistutorials.com/en/docs/3/handling_invalid_geometries.html
I followed all the points up to here:

As the algorithm runs, you will see a warning displayed in the Log
tab. 1 feature in the input layer had invalid geometry and was skipped
during the processing. The default setting to handle invalid geometry
in the Processing Toolbox is located at Settings ‣ Options ‣
Processing ‣ General ‣ Invalid features filtering and is set to Skip
(ignore) features with invalid geometries. This is a good default
setting, but if your input is large, you may miss this warning and may
not know that an input feature was skipped. You may want to change the
value to Stop algorithm execution when a geometry is invalid.

The next step describes what should happen

Back in the main QGIS window, you will see a new layer Colored added
to the Layers panel. Notice that the new layer is missing a state
which had invalid geometry. We now know that this particular state
polygon had invalid geometry but we don’t know what was the cause. We
can easily find that out. Search for and locate the Vector geometry ‣
Check validity algorithm.

I am getting errors as they described, here is a screenshot from what I'm doing:

But the layer does not get added to the list of layers. I'm using QGIS 3.18.0-Zurich (64bit).
For comparison I also ran Processing Tools -> Vector geometry -> Check validity, and that, after a long while, created the layer with highlighted places where geometry is flawed.


texts: Applying color ramp/map to mosaic using ArcGIS for Desktop?

texts:
My goal is to create a bold, colorful hillshade following this tutorial from the ESRI website.
Everything works great until I try to apply a color ramp or colormap. It somehow crunches the values and all I receive is this psychedelic hillshade.

I tried identifying my issue with the identify tool (how original) and there I noticed that no matter where I click on my mosaic, I receive the same output.

I must have botched up my values but don't know, when that happened. Maybe it is worth noting that I pretreated the DEM using focal statistics with a mean stat type.

EDIT:

Adding image of symbology.


texts: Installing Aggregate Polygons Function in PostgreSQL 9.3

texts:
I'm new to PostgreSQL and I've been trying to run a sql function written by Dr. Horst Duester, found in http://www.kappasys.ch/cms/index.php?id=65&L=5, which is used to aggregate detached polygons. It requires that two functions be installed first:

$ psql <dbname> -U <user> -h <server> -f cleanGeometry.sql

$ psql <dbname> -U <user> -h <server> -f aggregatePolygons.sql

both of which can be downloaded from the aforementioned web page.

I have managed to install the first, but when I try the second I get the following error message in the command prompt:

C:\Program Files\PostgreSQL\9.3\bin>psql -d test_db_1 -U postgres -h localhost -
f aggregatePolygons.sql

Password for user postgres:

psql:aggregatePolygons.sql:56: ERROR:  syntax error at or near "ο?Ώ"
LINE 1: ο?Ώ      ^

CREATE FUNCTION

psql:aggregatePolygons.sql:89: ERROR:  function _aggregatepolygonsfunction(geome
try, geometry, double precision, boolean) does not exist

GRANT

psql:aggregatePolygons.sql:93: ERROR:  aggregate aggregatepolygons(geometry, dou
ble precision, boolean) does not exist
psql:aggregatePolygons.sql:100: ERROR:  function _aggregatepolygonsfunction(geom
etry, geometry, double precision, boolean) does not exist

Has anyone encountered the same problem and found a solution?


texts: What is the naming convention used for US blockgroups codes?

texts:
I was looking at the naming convention for US blockgroups but could not find and reliable documentation for it. US Census website mentioned somewhere that the 5 digit US county code can be split as XXYYY where XX denotes the state code and YYY denotes an the county id within that state. Is there a similar convention for US blockgroups as well? I believe they use the FIPS system but wanted to get some more information on that.

Alternately, is there a source where I can get naming conventions for all these geospatial types?


texts: Convert Japan Administrative Boundary Data to Shapefile

texts:
I've downloaded Administrative Boundary data from the Japan MLIT (Ministry of Land, Infrastructure and Transport), but it is in an XML based format referred to as JPGIS or KSJ.

Are there any scripts, GDAL plugins, etc available for converting this format to a Shapefile or GeoJSON?


texts: Grouping layers in Table of Contents in ArcGIS for Desktop?

texts: I'm trying to group layers in the table of content so that I can expand or collapse them as needed.  Not sure how.  I'm using ArcGIS 10.3 for Desktop.


texts:
If I know that 10,000+ gallons of water were going to be dropped at a single point how could I model where it would flow.

I would prefer something that worked in ArcGIS and could create a layer.


texts: Reprojecting all geometries in PostGIS table?

texts: I created a spatial table with SRID:4326. Now I want to change total projection to SRID:32644 into a new table. The old table should remain unchanged.


texts: Customising a marker after performing a spatial search query

texts:
I'm have a geodjango app that uses leaflet for a front-end. I have integrated the search function on the front-end to be able to search a customer location in real-time, of which the code below does that perfectly and I'm happy. But the challenge now is after it displays the marker on the map, I'm unable to use the following code(see below) to customise the marker to show a pop up message ontop of the marker saying "I'm here" for example. i'm not sure what am I missing ?
<!-- search function for customer location -->
var searchControl = L.esri.Geocoding.geosearch().addTo(map);

var results = L.layerGroup().addTo(map);

searchControl.on('results', function (data) {
    results.clearLayers();
    for (var i = data.results.length - 1; i >= 0; i--) {
        results.addLayer(L.marker(data.results[i].latlng));
    }

<!-- Trying to customise the marker after the marker is placed on the map after performing a spatial search query above-->
L.marker(result).addTo(map)
      .bindPopup('Hi, I am here.<br> Which services are nearby?')
      .openPopup();


texts: ESRI-Leaflet updateFeature with a changed attribute or property

texts:
I'm updating a feature in ESR-Leaflet from a feature layer based on a call to a function. What I'm having trouble with is submitting this change to the feature service.

Currently I:
Get the feature based on it's ObjectID;
Reach into the feature.properties and change "SystemStatus";
Attempt to submit updateFeatures (and fail).

Code, and error below.

function ChangeStatus(objID){
    var retireFeatureID = hotspots.getFeature(objID);
    retireFeatureID.feature.properties["SystemStatus"] = "Deleted";
    hotspots.updateFeature(retireFeatureID, function (err, response) {
        if (err) {
             console.log(err);
             return;
        }
            console.log(response);
    });
    //hotspots.updateFeature(retireFeatureID);
    console.log("Updated the status of feature: " + objID);
}

The error I'm receiving is:

  {code: 400, message: "Cannot perform operation. Invalid operation
  parameters.", details: Array(2)} code: 400 details: (2) ["'features'
  parameter is invalid", "Object reference not set to an instance of an
  object."] message: "Cannot perform operation. Invalid operation
  parameters."

I'm not sure what could be causing my issue here.
Suggestions welcome.


texts:
I have been following a post here to create vector tiles direct from my database. My data is in WebMercator.

I am using the query:

SELECT ST_AsMVT(tile) FROM
    (SELECT dh_id,
        ST_AsMVTGeom(geom_go_mer,
           ST_MakeEnvelope(%(xmin)s,%(ymin)s,%(xmax)s,%(ymax)s, 3857), 4036, 256, true)
                AS geom FROM temp.buildings_google_mercator
                WHERE geom_go_mer && (SELECT ST_MakeEnvelope(%(xmin)s,%(ymin)s,%(xmax)s,%(ymax)s, 3857))
                ) AS tile

in Postgis 2.4.4 and parametising it with the webmercator coordinates calculated from the z,y,z tileschema (see here).

I get results, but as I zoom in, there are visibly more boundary disjoints. I looked at the resulting tiles, and it seems that the tiles are not calculated as I would expect.

In this diagram, the pink is the original geometry, the purple and green are two different adjacent tiles. Not only is the resulting tile geometry different from the original, but adjacent tiles do not overlap at all - the line within the green polygon should overlap the line within the purple polygon below.

I am confused whether

this is expected bahaviour for MVT
a bug / feature in Postgis MVT, or
something wrong with my interpretation of it.

I haven't seen bahaviour like this for other MVT maps, so I am not sure it's no. 1.


texts: How to check spotheight in closed polyline contour data using ModelBuilder or ArcPy?

texts:
I'm new in using ModelBuilder or Python script. So, i have a set of contour data and a good contour data should contain a spotheight* in every closed polyline contour data.

Is there any solution to check that every closed polyline of contour data should contain spotheight which also have same elevation or above.

*) spotheight is a point that located in closed polyline contour data which show the elevation value


texts: Interpolate Multiple Data With Kriging In ModelBuilder?

texts:
I am having problems building a model for running interpolation on multiple data using the Kriging method.

I have followed these steps:

Convert multiple data to geodatabase file
Used iteration on these feature classes
Plot Kriging tools in model

One of the inputs for the Kriging tool is a Z Value, it does not appear in the options, how can I include it?

I have attached my model

I am running in ArcMap 10.2

i just create new model concept and still Z value is empty


texts: Make Query Layer - field not available as unique identifier when using NVL

texts:
I need to create a spatial view with Make Query Layer. My query includes outer joins, so I must combine every table's ID field to create a unique identifier. Since there might not always be a child for every parent record, I use NVL to replace null occurences in the ID fields with a constant value (say 9999).

My problem is that these fields using NVL are not visible in the Unique Identifier Field(s) (optional) list parameter.  Any idea why?

Or is there any other method I can use to get unique IDs for my query layer?

I am not allowed to create my view directly in the database.

I use ArcGIS 10.3.1.  and the underlying database is Oracle 12c. All ID fields are integers (NUMBER(10, 0)).


texts: Trying to recode raster with CSV using the Grass r.recode.attr but getting error "name 'xrange' is not defined"

texts:
I'm trying to recode a large GSSURGO soil raster with CSV using the Grass r.recode.attr following these instructions https://pvanb.wordpress.com/2014/12/13/recode-your-raster-file-in-grass-gis-using-a-csv-file/ but I'm getting error "name 'xrange' is not defined"

When I looked it up I found that this issue may be because "xrange" changed to "range" in the transition from from Python 2 to Python 3. The instructions for the r.recode.attr module are from 2014 and I can't find much else out there the talks about r.recode.attr. However, it seems like a very valuable tool.

My raster, instead of being classed by values has a unique key which only has meaning if you can reclassify it from a related CSV table. Does anyone know of any alternative using GRASS or QGIS?


texts: Is there GIS Video "that Mom can understand"?

texts: Working as a professional within the GIS discipline, I often find it challenging to describe to people what GIS is (in a face to face conversation).  Telling them that GIS stands for Geographic Information Systems, which encompasses a set of technologies that can be used to analyze, map, and organize spatial data is a bit over most people heads.  On the contrary, telling people that GIS can be used to make maps on a computer falls far short of its many potential uses.  With this, I'm looking for some effective video examples that one could show somebody that explain both in words and visually what GIS is, and how the technology is utilized (something your Mom could understand :))?


texts: Is Google Earth Engine SNIC Segmentation algorithm inconsistent?

texts:
I have been using GEE for about a month now and I ran into some issue

I am trying to use OBIA classification based on this presentation by Sai Cheemalapati, EEAI 2019.

The main problem is when I perform SNIC algorithm on sattlite image, everytime I change scale, it recalculate and give me a new result. I tested my hypothesis that SNIC is sensitive to scale by exporting region in 2 scales. And to my surprise, both image did show inconsistency (large area segmentation on scale of 10 vs small area segmentation on scale of 2)

Is it possible to specify scale to make SNIC algorithm segmentation on that scale only?

Is this a bug?

Or I am missing something in the documentation, and if so can you advise?

Example Image

SNIC with scale size 2

SNIC with scale size 10

The simplified demo code can be found here

var S2 = ee.ImageCollection("COPERNICUS/S2_SR"),
    ROI = ee.Geometry.Polygon(
        [[[10.093076096894833, 55.45897533022183],
          [10.093076096894833, 55.407161055045385],
          [10.231778489472958, 55.407161055045385],
          [10.231778489472958, 55.45897533022183]]], null, false);

// Prepare Image
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

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

var image = S2
  .filterDate("2019-03-01", "2019-04-01")
  .filterBounds(ROI)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 80))
  .map(maskS2clouds)
  .median()
  .clip(ROI)

// Prepare SNIC
var kernel = ee.Kernel.gaussian(3);
image = image.convolve(kernel);

var seeds = ee.Algorithms.Image.Segmentation.seedGrid(36);
var snic = ee.Algorithms.Image.Segmentation.SNIC({
  image: image,
  compactness: 0,
  connectivity: 8,
  neighborhoodSize: 64,
  size: 3,
  seeds: seeds
});

// Visualize
Map.addLayer(image, {
  bands: ["B4", "B3", "B2"],
  min: 0,
  max: 0.255
}, "RGB")
// Map.addLayer(snic.randomVisualizer(), null, "snic")

// // Run this on 1st time
// Export
Export.image.toAsset({
  image: snic.select(["clusters"]),
  description: "cluster_scale-10",
  assetId: "cluster_scale-10",
  region: ROI,
  scale: 10,
  maxPixels: 1e13
})
Export.image.toAsset({
  image: snic.select(["clusters"]),
  description: "cluster_scale-2",
  assetId: "cluster_scale-2",
  region: ROI,
  scale: 2,
  maxPixels: 1e13
})

// // For other time
// var user = "your-user-name"
// var snic_s_2 = ee.Image("users/" + user + "/cluster_scale-2");
// var snic_s_10 = ee.Image("users/" + user + "/cluster_scale-10");

// Map.addLayer(snic_s_2.randomVisualizer(), null, "SNIC 2 m")
// Map.addLayer(snic_s_10.randomVisualizer(), null, "SNIC 10 m")


texts: Why does Raster Calculator return GRID ERROR-Unable to open input grid?

texts:
I'm calculating Simple Ratio (SR) from NDVI using Raster calculator in Arcgis. Here is expression I wrote in Raster calculator.

SR=(1+float[NDVI_apr.tif])/(1-float[NDVI_apr.tif])

But when I evaluate this expression, the following error always appear on the screen.

GRID ERROR-Unable to open the input grid: 1+FLOATC:\document\...
Error in executing grid expression.

What is wrong with this expression?


texts: how to measure polygon perimeter / Distance in openlayers 3?

texts:
i am using openlayers 3 measurement tool for measuring line lengths and polygons areas, and this example worked for me: http://www.rhone-mediterranee.eaufrance.fr/milieux-aquatiques/poissons/js/openlayers-v3.19.1/examples/measure.html ,
but i want to add something like getlenght or get perimeter for the polygon so when i draw a polygon i get the area and the perimter also ??


texts: What is the efficient way to dissolve polygons using Spatialite 3.0

texts:
I have polygon table called "tableA" with 72000 records. It has a geometry column called "geometry", and a text column called "field1". I want to  aggregate\dissolve all polygons with the same "field1" value.

I have a Spatial Index in Geometry and an index in field1.

I tried this:

Select f.field1 as field1, st_union(f.geometry) as geometry
From tableA as f
Group by field1;

and its taking too long, I had to cancel it after being processing for 1 hour. Using Arcgis it toke me 5 minutes, so I must be doing something wrong.

So, is there a better way to preform this operation using spatialite? Is the Spatial Index being used this way?


texts: finishCondition based on feature being drawn

texts:
I have a draw interaction in Openlayers 4.

I'm using it to draw polygons. I would like to continue drawing the polygon if its area is less than a default value.

I tried using the property finishCondition, but it takes an ol.MapBrowserEvent and feature is not accesible.

finishCondition(event) { /* return (true or false) based on feature */ }

It seems that the event is fired before feature is added to the source so I can't take the feature from the source.

Is there any way to do this?