texts: How to turn off Google map data copyrighted window in Google Street layer (Qgis - Openlayers Plugin)?

texts:
Possible Duplicate:
  How to close watermarks in openlayers plugin

After a few weeks not using Qgis, this morning I opened it up, accepted a request to update Openlayers Plugin. Right now, I found a window saying "Map data ©2012 Google, Mapabc" at the center of the screen, following me every where I drag the map. Any one experienced this, please tell me how to remove that screen?

I am using Qgis 1.7.4 on Windows XP, Openlayers Plugin version 0.90.


texts:
SRTM3 .hgt files have 1201*1201 DEM elevation data, I extract .hgt files data elevation into file Elevation.json with index.js with the help from this GitHub project.

Then I rename Elevation.json into Elevation.txt, replacing all ],[ from inside the converted .hgt file in json array into ]^l[ with MS Word to get each row from json at new line.

Finally, I import this .txt file into .xlsx file with "comma" as separating delimiter.

I expect as final result in Excel, just 2101 rows and 2101 columns.

Why is there an error always at these cells?

I shortcuted the first and the last 3 columns and the first and the last 2 rows to show you the result

.

Note: blank cells = no found data.

I tried different .hgt file from world wide, the problem still exist.


texts: Editing ArcGIS metadata elements using Python?

texts:
I have a shapefile with an xml metadata file. Most of the elements are pre-populated using the previous editions xml file. However, I am trying to automatically update the current edition within a larger script.

I cannot update my Python library with any new modules (for example the metadata module seen in other questions) as I am using a work computer which is locked down and cannot install anything.

The three elements that need updating are the title, the revision date and the edition date.

I have managed to update these using the tree and an index to find the tags.

However, how do I search for the metadata element using the tag title? So for example, the the edition date tag is "redEdDate".

In addition, will the metadata tree index ever change if/when additional information is added? Or do the trees have a set format, therefore the same indexes. Im worried that any update to ArcGIS metadata, will affect the tree, and therefore the index for the specific tags mentioned above, and cause error within the script.

The issue above is shown via the script below. Both datasets use the ArcGIS metadata format. However the same tags have different indexes. For example, the MCMS Polygon dataset uses an index of "root[4][0][0]" for the edition date. However the MCMS Exclusion Zone edition date has an index of "root 3[0][0]".

My script so far is:

import arcpy, os, sys, datetime
import xml
import xml.etree.ElementTree as ET
ws = arcpy.env.workspace = r"path/to/folder"
today = datetime.date.today()
date = today.strftime("%Y%m%d")

#Update the MCMS polygon metadata titles, update date and edition date
for f in os.listdir(ws):
#Find the polygon xml file
    if f.endswith("A.shp.xml"):
        fpath = os.path.join(ws, f)
#Identify the metadata tree
        tree = ET.parse(fpath)
        root = tree.getroot()
#Set the title and date variables to the relevant metadata tag index
        editiondate = root[4][0][0]
        reviseddate = root[4][0][5][0]
        title = root[4][0][7]
#Update the tags with the new data
        editiondate.text = today.strftime("%Y-%m-%d")
        reviseddate.text = today.strftime("%Y-%m-%d") + "T00:00:00"
        title.text = "MCMS (polygon)"
#Write the updates to the xml file
        tree.write(fpath)

#Update the MCMS exclusion metadata titles, update date and edition date
for f in os.listdir(ws):
#Find the exclusion zone xml file
    if f.endswith("Zones.shp.xml"):
        fpath = os.path.join(ws, f)
#Identify the metadata tree
        tree = ET.parse(fpath)
        root = tree.getroot()
#Set the title and date variables to the relevant metadata tag index
        editiondate = root[3][0][0]
        reviseddate = root[3][0][5][0]
        title = root[3][0][7]
#Update the tags with the new data
        editiondate.text = today.strftime("%Y-%m-%d")
        reviseddate.text = today.strftime("%Y-%m-%d") + "T00:00:00"
        title.text = "MCMS Exclusion Zones"
#Write the updates to the xml file
        tree.write(fpath)


texts: How do I get the attributes of a particular polygon in qgis

texts: I downloaded a shapefile that has demographic data for census tracts for a city -- each represented by a polygon in the vector layer. I can click the attribute table icon to see the whole sheet of attributes for each polygon -- but is there any way to click a particular polygon and see it's attributes. I tried right clicking and double clicking.


texts:
I am a QGIS new user and I want to extract the buildings height of a territory. I clipped the raster with the LiDar informations with the buildings polygons; after I have extract the contours level inside each polygon.

Now I want to extract the maximum value (top of the roof) and the minimum value (base of the roof) and assign them to each polygon as attribute. How can I do it?


texts: Can't get array_foreach function working

texts:
I have a STRING field with values. When I want to apply some additional text after each value using an array function, I get the text added to the end of the whole initial string instead of each value.
I tried space separated values, coma separated values, coma separated values with spaces but nothing is working.
Values can be numbers or letters, like 1 5 XY for example.
I applied the following expression in a virtual field calculator:
array_to_string( array_foreach( array("values"),@element || 'sometext' ), ' ' )

QGIS returns me:
1 5 XYsometext

instead of:
1sometext 5sometext XYsometext


texts: QGIS Add Vector Geometry to Polygon query

texts: I have intersected a landcover map with 10km grid squares and then used the vector geometry tool (Add/ Export geometry) to calculate the area of each landcover feature within the grid square. In the resulting attribute table I have four new columns - area, percent, area2 and perimeter. Could someone please tell me what the area2 value represents? I have had a look online but can't find any answers.


texts: How to produce maps with copyright imagery on it in GRASS?

texts:
I'm having some trouble with a GIS that I've got to produce.

I work for a contractor and we'll have to produce maps to do the work for our client. The client has specified that if we print maps out (which we will have to) they will have to have various images on - the client's logo, map provider's logo, copyright information etc. We do have maps provided for us in PDF which have these logos on but they are an absolute pain to use, zoom in to a decent scale and keep the logos on. Using a GIS like GRASS would be much easier.

I'm new to GRASS and haven't found the option to add images to printouts yet. I've been messing around with the Cartographic Composer but haven't managed to find the option to add images - just text or automatically-generated content like legends.

Can anyone tell me where to look, or whether later versions of GRASS have this?


texts: Is it possible to paste features into QGIS from another application?

texts:
One of the coolest features in my QGIS workflow is the ability to copy features (select and copy) and then paste the results into, say, Excel or Openoffice - you get your attributes table plus a WKT column that describes the geometry:

wkt_geom               id    name
POINT(250000 9012000)  1     Home

(You can then parse this directly in Excel using something like =MID(A2, 7, 6) to get the X value of the point).

What would be even better is to be able to edit that table in the external app, copy it again, and paste it back into QGIS (with the relevant layer in edit mode). However, I've tried that without success - QGIS remembers its set of copied features from the previous operation, ignoring the system clipboard (which I've confirmed has the right table stored).

Any clever ways to make this work?


texts: Checking whether two rasters intersect using R?

texts:
I would like a method for checking whether two rasters (or a raster and an extent) overlap that returns a Boolean and doesn't involve an if/then statement. Something like the gIntersects function. Here's my current approach:

#load libraries
library('raster')
library('rgeos')

#extent rectangle
rect <- extent(956328.8, 1036613, 825238.9, 936911.5)

# build blank raster
r <- raster(xmn=1035792, xmx= 1116792, ymn=825303.6, ymx=937803.6, resolution = 100)
r <- setValues(r, 0)

#check whether raster overlaps extent
r.crop <- crop(r, extent(rect))
if(is.null(r.crop)) {
    print('no overlap')
}else{
    print('overlap')
}


texts: arcpy and multiprocessing: switch workflow to shapefiles to avoid schemalocks

texts:
I need to processe a large dataset  through a bunch of interwoven scripts using arcpy. On a small test dataset my workflow functions fine so far but when it comes to the large dataset i see no other way as to create subsets to iterate on like below:

jobs = []
for i in list_of_subsetFCs:
    args = (i[0],i[1],10,5,1,"string") #inputFC, parameters and outputfile syntax
    p = multiprocessing.Process(target=my_function, args=args) #my function is importet
    jobs.append(p)
    p.start()

my_function makes use of:

FieldMappings
MakeFeatureLayer_management
GetCount_management
Select_analysis
Intersect_analysis
SplitLineAtPoint_management
SpatialJoin_analyses
AddField_management
Append_management

Is it possible and/or desirable to change the whole workflow in my_function to shapefiles in order to avoid GDB schema-locks or is there perhaps another possibility?

(Arcgis 10.2.2, 64bit background processing)


texts: Creating drive time polygon around map point?

texts:
I'd like to create a 30 min drive time polygon around 400 separate points I have plotted to google maps.

I have discovered http://www.freemaptools.com/how-far-can-i-travel.htm but this seems to freak when you get near the sea (many of my points are).

And also http://cartoo.dyndns.org/ which works well, but I can't then extract the polygons in .KML form.

Anyone come across any tools able to help me?


texts: When calcuating distance between points on Earth why are my Haversine vs. Geodesic calculations wildy diverging?

texts:
I am getting wildly diverging distances using two approximations to calculate distance between points on Earth's surface. I am using the Haversine (vectorized) approximation (spherical earth) and the more precise (presumably) geopy.distance.geodesic(ellipsoidal earth) .  See here for discussion of Haversine n Geodesic

As you can see I am off by five percent as the distances between points becomes large. Is this divergence due to rounding error in Haversine? Do I indeed trust the Geodesic? Here is code:

import numpy as np
lat = np.linspace(35,45,100)
lon = np.linspace(-120,-110,100)

data = pd.DataFrame({'Latitude':lat,'Longitude':lon})

def Haversine(v):
    """
    distance between two lat,lon coordinates
    using the Haversine formula. Assumes one
    radius. r = 3,950 to 3,963 mi
    """
    from timeit import default_timer as timer
    start = timer()
    R = 3958 # radius at 40 deg 750 m elev
    v = np.radians(v)

    dlat = v[:, 0, np.newaxis] - v[:, 0]
    dlon = v[:, 1, np.newaxis] - v[:, 1]
    c = np.cos(v[:,0,None]
    a = np.sin(dlat / 2.0) ** 2 + c @ c.T * np.sin(dlon / 2.0) ** 2

    c = 2 * np.arcsin(np.sqrt(a))
    result = R * c
    print(round((timer() - start),3))
    return result

def slowdistancematrix(data):

    from geopy.distance import geodesic
    distance = np.zeros((data.shape[0],data.shape[0]))
    for i in range(data.shape[0]):

        lat_lon_i = data.Latitude.iloc[i],data.Longitude.iloc[i]

        for j in range(i):

            lat_lon_j = data.Latitude.iloc[j],data.Longitude.iloc[j]

            distance[i,j] = geodesic(lat_lon_i, lat_lon_j).miles
            distance[j,i] = distance[i,j] # make use of symmetry

    return distance

distanceG = slowdistancematrix(data)
distanceH = Haversine(data.values)

plt.scatter(distanceH.ravel(),distanceG.ravel()/distanceH.ravel(),s=.5)
plt.ylabel('Geodesic/Haversine')
plt.xlabel('Haversine distance (miles)')
plt.title('all points in distance matrix')

I would rather use the vectorized version becuase it is fast. However,the 5% is too big for me to be comfortable with it. Supposedly Haversine is only suppose to be off by .5%.

Update
Solved. Matrix algebra error.


texts: Tools to help migration from Mapserver to QGIS Server

texts:
I'm going to migrate a Mapserver WMS server to QGIS Server.

I'd like to know if there are some tools or plugins that convert Mapserver mapfile (.map) to QGIS project file (.qgs). WMS responses (GetCapabilities, GetStyles...) to QGIS project file would be also be a correct solution.

Ideally it would take care of :

WMS server metadata
layer's name
layer's title
layer's style
layer's HTML template
layer's minscale/maxscale view
layer's groups
layer's metadata
....


texts: Simple field calculator equation in ArcGIS 10 not working

texts:
I'm new to ArcGIS 10 and am trying to do a really simple field calculator equation and it is not working. I'm just trying to find the percentage of a certain demographic trait by doing:

[asian]/[total]*100

It is populating about 10% of the column correctly and then I get an error:

  There was a failure during processing, check the Geoprocessing Results
  window for details

Obviously my syntax is not incorrect or else it wouldn't start populating the field correctly.

Does anyone have any idea why this is happening?


texts: ArcGIS FeatureService/Dynamic Map Service using SQL joins and aggregations

texts:
I have a bunch of operational tables, some contain spatial data, some do not. In order to fulfill all my operational mapping reports, I will need to have the flexibility to have feature services or dynamic map services that perform table joins and various aggregations on the underlying geodatabase (or client side). All the tables are in the same enterprise geodatabase instance but not all of them are spatial tables (include a geometry column).

What are my various options to do this? So far I have come up with:

Joining and Relating tables (http://help.arcgis.com/en/arcgisdesk...0000002n000000)
It's my understanding, however, that feature services do not support data from joined tables
(http://resources.arcgis.com/en/help/...00000007000000)
Query layers
Client side Feature Service joins and aggregations
Custom Geoprocessing service (ArcPy) that returns a FeatureSet made up of a custom FeatureClass (input schema?)
generated by running the join,aggregation query server-side, creating the custom FeatureClass list and returning that to
the client, the idea being that I would make the custom geoprocessing service look like a FeatureService or a Dynamic
Map Service to the client.

I'd like to get people's opinions in terms of maximum flexibility and scalability ( I don't want to run into a situation where the client must handle large numbers of records in order to perform client side aggregations and joins, that stuff should be done up front in the server and the aggregated data should be made available to the client. )


texts:
Is there a method to save GDAL/OGR Utilities --config default options?

For example I almost always want to use the gdaladdo option for smallest possible jpeg-in-tiff overviews, but that's a lot of frequent typing and/or pasting:

--config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL

Can I set my shell to always use those options unless I specifically want otherwise? (Windows CMD most of the time, but Linux answer would be good too).


texts: Raster calculator giving a different result in processing modeler Qgis 3.2.3

texts:
I'm running Qgis 3 processing modeler and I am trying to invert the values of a  raster file.

I had previously asked this question and the answer I came up with works with the raster calculator when processing the layer.  However when I try the same thing within a processing modeler I get a blank raster file.

This is the raster info from the proximity output.

Driver: GTiff/GeoTIFF
Files: C:/Users/Timothy/AppData/Local/Temp/processing_54c8e3d84eb144a88f3f486c6ae75076/c5fc54eb536f4b1cbe57b7ef55506bb2/gdal_proximity_1_proximity.tif
Size is 1000, 1000
Coordinate System is:
PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",37.5],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-2249531.050643497612327,289551.660125185269862)
Pixel Size = (0.577937172434758,-0.330446977802843)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-2249531.051,  289551.660) (121d56'57.53"W, 37d19'33.78"N)
Lower Left  (-2249531.051,  289221.213) (121d56'53.88"W, 37d19'23.56"N)
Upper Right (-2248953.113,  289551.660) (121d56'34.71"W, 37d19'38.79"N)
Lower Right (-2248953.113,  289221.213) (121d56'31.06"W, 37d19'28.56"N)
Center      (-2249242.082,  289386.437) (121d56'44.30"W, 37d19'31.17"N)
Band 1 Block=1000x2 Type=Float32, ColorInterp=Gray


texts: Mapserver - Kernel Density Layer - weight

texts:
I have created a map file broadly as shown below. I then created 6 points as in the DATA_TABLE. The 3 on the left have val 100 and the three on the right have val 1.

When I view this layer it renders as follows. With both sides looking the same? Why is the val value not affecting the weight and hence the colour ramp.

Have I misinterpreted something?

I was expecting the higher valued items on the left to be "redder".

LAYER
  NAME feature_heatmap
  TYPE raster
  CONNECTIONTYPE kerneldensity
  CONNECTION "point_data"
  STATUS ON

  PROCESSING "RANGE_COLORSPACE=HSL"
  PROCESSING "KERNELDENSITY_RADIUS=50"
  PROCESSING "KERNELDENSITY_COMPUTE_BORDERS=ON"
  PROCESSING "KERNELDENSITY_NORMALIZATION=1"
  OFFSITE 0 0 0

    CLASS
          STYLE
            COLORRANGE  "#0000FFFF"  "#FF0000FF"
            DATARANGE 0 255
          END
    END

 END

 LAYER
      PROJECTION
        "init=epsg:4326"
        "init=epsg:3857"
        "init=epsg:27700"
      END # PROJECTION

      NAME "point_data"
      STATUS ON
      TYPE POINT
      CONNECTIONTYPE POSTGIS
      <OBFUSCATED CONNECTION DETAILS>
      DATA "wkb_geometry from (SELECT wkb_geometry, feature_id as id, val::INTEGER as VAL FROM DATA_TABLE)))  AS FOO USING UNIQUE id using srid=4326"
      CLASSITEM val
          CLASS
            STYLE
              SIZE [val]
            END
          END

 END


texts: Arcpy doesn't recognize spatialite tables after editing with QGIS

texts:
I'm working on a project that would like to use spatialite as a go-between for field edits using QGIS/IntramapsRoam and then loading the data to a file geodatabase.

The workflow is as follows
Spatialite (created in ArcGIS using 'SPATIALITE' geometry) -> QGIS -> ArcGIS 10.2

I'm running in to a compatibility issue with arcpy and spatialite that renders my workflow impossible and am wondering if this has been encountered before and if there are any known solutions/workarounds.

#Create the database
sqllite_database = r'.\test.sqlite'
arcpy.gp.CreateSQLiteDatabase(sqlite_database, “SPATIALITE”)
arcpy.CreateFeatureclass_management(r'.\test.sqlite',
                                    'points', "POINT")
arcpy.AddField_management(r'.test.sqlite\main.points',
                          "num_field", "TEXT")

After creating the database, I can access it with arcpy.

I then add a point in QGIS and save edits.

I can no longer access the layer with arcpy.

Oddly enough, I can still access the table from ArcMap, but arcpy is the functionality that I really need.

Given the lack of informative error messages and documentation for arcpy with spatialite, does anyone know what might be causing this? Any suggestions for solutions or workarounds?

Update: Here is a screen showing the layer drawing in ArcMap and returning False in arcpy. This seems like errant behavior.

Update 2: I tested again using ArcGIS10.3.1 and it still failed. Apparently you cannot edit sqlite databases in QGIS and retain arcpy functionality.


texts: PostGIS ST_Within ST_intersects differing SRIDs

texts:
I have a layer of points and want to join some data from polygons where the two coincide. I might have several points in one poly.

The points were generated from a CSV as;

UPDATE t_osab_with_voa_ndr SET geom_point_4258 = ST_SetSrid(ST_MakePoint(longitude, latitude),4258);

Where they are stated as being ETRS89 projection which is SRID 4258

I'm running into problems as the points have SRID 4258 and the Polys are 27700 I've tried using ;

ALTER TABLE t_osab_with_voa_ndr
 ALTER COLUMN geom_point_4258 TYPE geometry(point,27700)
  USING ST_Transform(geom_point_4258,27700);

When I look in QGIS I see this;

But when I run;

create  or replace view v_spatialjoin_northeast_st_within1 as
select * from t_osab_with_voa_ndr , d_geom_northeast  where st_within(t_osab_with_voa_ndr.geom_point, d_geom_northeast.geom);

I get no results when I can see I should

Question 1.
Can i combine the ST_Within / ST_Interects and specify the SRID in the process or am i better transforming the points first?

EDIT:
Reason for my blanks was SQL problem with ST_Within, but still interested to know if you can set SRID or what is the fastest way to execute this


texts: Export lon/lat from vector. GRASS GIS / QGIS

texts: I have some data as .shp (polygone representing Antarctica continent). I would like to know if there is a way to create a raster containing the lat/lon for each grid point ... Has somebody here tried to to that?


texts: Using ArcGIS Polygon to Raster tool does not divide resulting raster on cells?

texts:
I try to create raster composed of small (5x5 m) cells (pixels) from a polygon in ArcGIS-10.
I created a polygon (as a single object), and use Feature to Raster Tool and Polygon to Raster tool (of Conversions tools). I assigned "output cell size" as 5 m.
In both tools the program creates raster but without cells, so as a single object. The attribute table of the final raster has 1 line with Value = 1 and  Count = 61394.

So, as I understand, it has 61394 cells with value = 1, but for some reason the program does not divide it into the cells themselves.

My first polygon and the data frame have UTM Coordinate system, so all data - in metres.

What I have to do for create this raster correctly?

My ultimate goal -  to cut my first big polygon (that was created as a single object) by grid of equal small cells 5x5 m. I went toward the creation  raster at first, because I didn’t find a direct way.


texts: Copy geometry data as HEXEWKB from csv to table in PostGIS

texts:
I have a table in PostGIS 2.3 (PostgreSQL 9.5, and pgadmin III) where I store Multipolygons. I have a csv which is a straight export from another database, including the geometry column in hexewkb format (as I understand). I cannot figure out a way to simply insert this data into the table. Here is an example :

CREATE TABLE tout.battest(
 id_battest serial NOT NULL,
 geom geometry(MultiPolygon,2154),
 CONSTRAINT prk_constraint_battest PRIMARY KEY (id_battest));

 INSERT INTO tout.battest(geom) VALUES ('010600000001000000010300000001000000050000007041F528CB332C413B509BE9710A594134371E05CC332C4111F40B87720A594147E56566CD332C4198DF5D7F720A594185EF3C8ACC332C41C03BEDE1710A59417041F528CB332C413B509BE9710A5941');

This insert into raises the error "Geometry SRID (0) does not match column SRID (2154)". So I tried this:

INSERT INTO tout.battest(geom) VALUES (st_setsrid('010600000001000000010300000001000000050000007041F528CB332C413B509BE9710A594134371E05CC332C4111F40B87720A594147E56566CD332C4198DF5D7F720A594185EF3C8ACC332C41C03BEDE1710A59417041F528CB332C413B509BE9710A5941',2154));

Which then raises an error like "the function st_setsrid(unknown, integer) is not unique"

I then tried to call various functions, but without any success, sometime even getting an error "function does not exist" when I actually find it in the function list of my public schema. Any idea on how to do what I want?


texts: Adding Feature IDs to existing MbTiles or on Mapbox

texts:
Question TLDR:
Is there a way to add feature IDs to an mbtiles file after it is generated?
Or alternatively: is there a way to add feature IDs to a source via mapbox either client side or otherwise?
Details:
I have a large geoJson dataset of over 20Gb that I exported to mbtiles with tippecanoe which took around 3 days in order to upload to Mapbox as a tileset.
Unfortunately I forgot to set the --generate-ids flag when exporting, so my features have no feature ids and therefore do not have the functionality I require on mapbox (setFeatureState on particular features).
Is it possible to add in featureIds after the fact to mbtiles files?  On the mapbox side I am aware of both generateId (unfortunately geojson specific) and promoteId here which I could use with my data but unfortunately does not work because the source is being created with a URL and gives the error:
sources.<source-name>.promoteId: a source with a "url" property may not include a "promoteId" property

Otherwise, I will have to re-generate the entire mbtiles file which seems like a bit of a waste.


texts: "Selected" button on the attributes window

texts: I'm following some old instructions that are not up to date with the current version of ArcMap. I'm working on a table that is a layer. When I open it the instructions state that I have to click the 'selected' button on a specific attributes window. What button are they referring to?


texts: Converting Python (ArcPy) script into tool or toolbox script?

texts:
I've written a simple python script that allows me to quickly create a copy of a template Mxd, pointing to different shapefiles. The user specifies the new location, a name for the new mxd, then the folder that the shapefiles are stored in and the name of a point, polyline, and polygon shapefile. This script works fine.

newFilePath = r"X:\folder\outputfolder"
mxdName = "New_MXD"
shpFolder = r"X:\folderpath\templates"
shape_points = "input_P"
shape_polylines = "input_L"
shape_polygons = "input_G"

template_points = r"X:\folderpath\template_P.shp"
template_polylines = r"X:\folderpath\template_L.shp"
template_polygons = r"X:\folderpath\template_G.shp"

import arcpy
MyMXD = r"X:\folder\template.mxd"
mxd = arcpy.mapping.MapDocument(MyMXD)
outMxd = newFilePath +"\\"+ mxdName+".mxd"
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
    if lyr.dataSource == template_points:
        print lyr.dataSource
        lyr.replaceDataSource(shpFolder, 'SHAPEFILE_WORKSPACE', shape_points)
        print "data replaced"
    if lyr.dataSource == template_polylines:
        lyr.replaceDataSource(shpFolder, 'SHAPEFILE_WORKSPACE', shape_polylines)
    if lyr.dataSource == template_polygons:
        print lyr.dataSource
        lyr.replaceDataSource(shpFolder, 'SHAPEFILE_WORKSPACE', shape_polygons)

mxd.saveACopy(outMxd)

I would like to turn this script into a tool that my colleagues can run without having to open the script. I've seen Esri guides for adding scripts to a toolbox, and for creating python toolboxes, but its not clear which I need to do. Any suggestings would be great.

I've tried the first option, and replaces the first block with:

NewFilePath = arcpy.GetParameterAsText(0)
mxdName = arcpy.GetParameterAsText(1)
shpFolder = arcpy.GetParameterAsText(2)
shape_points = arcpy.GetParameterAsText(3)
shape_polylines = arcpy.GetParameterAsText(4)
shape_polygons = arcpy.GetParameterAsText(5)

When I tick "Run Python script in process" the script runs but nothing happens. When I untick this and tick "show command window when excecuting script" the script fails and I get ERROR 000714: Error in script ViewSismageExport.

Is setting up parameters more complicated than this? I first tried setting the folders to be folders and the shapefiles to be shapefiles, but I've not got everything as string, with it still not working.

Does anyone know what I'm doing wrong?


texts: Floating point to 8-bit unsigned raster conversion using ArcGIS for Desktop or ERDAS IMAGINE?

texts:
I have an image with value ranging from -1 to +1.

I want to convert this into 0-255 (8bit).

How can this be done using ArcGIS for Desktop or ERDAS IMAGINE?


texts: Position of vehicles relative to lead vehicle

texts:
I have a set of vehicles (think railway) that are in a (semi)fixed formation. I know the lead vehicle and I know the GPS coordinates of the lead vehicle at any given time. Sometimes, I can rely on the on board systems on the vehicle to report the other vehicles in the set, along with the order of the vehicle. This information is not always available though, depending on the type of vehicle.

I have already got a pretty good method of determining which vehicles are travelling together. What I need to solve is in what order are they travelling, starting with the lead vehicle and working backwards.

I have limited spatial data to work with, but what I do have is vehicle number, timestamp, latitude and longitude. The data is stored in MySQL, the coordinates are stored both separate and as a Point(long, lat) for spatial indexing.

I think I'm getting close, around 80% accuracy, based on a small sample. What I am basically doing is measuring the distance from the known lead vehicle at any given timestamp to the coordinates from the other vehicles using ST_Distance_Sphere. It's pretty good, but not quite there yet.

Has anyone done anything similar before? Or does anyone have any good ideas how I can tackle this in a better way? If possible, I'd prefer to use the GIS functions in MySQL, but this isn't a prerequisite.

Here is the map of the combined GPS points for three vehicles:

This is a simple diagram of the real life situation, with vehicle 1 being accurate, 2 and 3 being known but the positions unknown:

Update 1: Clarification

I have GPS coordinates and timestamps for all vehicles. What I need is the order in which the vehicles are travelling relative to the lead vehicle.

The SQL I already have gives me the definitive lead vehicle (based on some other data) and the other vehicles that are travelling with the lead. What it doesn't give me is the ordering of the vehicles (except the lead, or position 1).


texts: ArcGIS map export is exporting some features but not others

texts:
I feel this is pretty simple but I haven't had much success searching for an answer and I'm also running out of time so I'm resorting to asking here.

I am trying to export a map with the exact same settings as I always do, but it is choosing not to export the more complicated features such as contours, water and our vegetation inventory while successfully exporting roads, annotation and all the layout components. The problem does not seem to be confined to a certain geometry type. I get the same result no matter the format or settings. Some of the data is on a local server while most of it is from a shared corporate server.

The main culprit is probably the fact that I'm using Data Driven Pages, but I get the same result when I turn them off.

Any ideas? It's not even attempting to export them as it completes in a few seconds while normally I would expect it to take longer. I unfortunately do not have any graphics programs to mess around with.


texts: Remove ESRI logo from map with Silverlight API v2.1

texts: Is there a way to remove the ESRI logo from the map display with the ESRI Silverlight API v2.1.  This was possible in the 1.x version, but the old method no longer works.


texts:
We are preparing a Tileserver to serve tiles with a time dimension that updates every day. I was wondering what could be the best approach for that or the ideas we had are a good way.

To illustrate the problem:

We have observations for every full hour over 10 hours. That would be then 10 timesteps if implemented with a timeslider. The observations change every day.

Our idea was to make 10 layers which are accessible over a different layername. The layername would reflect the time of the observation(e.g. ISO 8601 timestring) or could be a simple number like 1, 2, 3, ... for every timestep with a corresponding meta.xml containing the date of the first layer and the interval between the timesteps.

Of course it would be the best if the meta.xml would be something like a WMS GetCapabilities. Different tile serving softwares as MapProxy or Mapcache serve some kind of metadata.xml for TMS and WMTS GetCapabilities for WMTS but seem unable to insert a time dimension. How could we implement this? There was a post asking for implementations of the WMS time dimension for a WMTS but maybe in the meantime there where some developments. There is also a blog post on how the guys at nasa added a Dimension block to the WMTS GetCapabilites but I don't know which software they used for it. In the end we would prefer to use some existing standards over implementing our own structure.

Is this a good way to go? Do you think there is a better solution?


texts:
I'm still pretty new to arcpy and am still learning. I'm working on a script that has a polyline as its input. I'm able to traverse the polyline and get points at certain intervals. These intervals are represented as segments If I have segment AB, where A and B are points on the polyline, how can I traverse the points on the polyline in between A and B and grab those points to add to a shapefile?

Thanks in advance


texts: Any efficient workaround to reshape and aggregate multi-layers raster grid in R?

texts:
I have historical daily minimum/maximum temperature observations data in two big netCDF file, so I used raster::stack to import those in RasterStack object. However, here I want to reshape and aggregate two different multi-layers raster grid into one plain text ASCII file (want to replicate of this).

reproducible data

Here is the reproducible example for historical daily minimum and daily maximum temperature observation in multi-layers raster grid:

T_max <- raster(xmn=5.75, xmx= 15, ymn = 47.25, ymx =55,res=c(0.25,0.25))
Deu_df_max <- do.call(stack,lapply(1:5479,function(i) setValues(T_max,round(runif(n = ncell(T_max),min = 2,max = 35)))))
names(Deu_df_max ) <- paste0('X',gsub('-','.',ymd('1980.01.01') + days(1:5479)))

T_min <- raster(xmn=5.75, xmx= 15, ymn = 47.25, ymx =55,res=c(0.25,0.25))
Deu_df_min <- do.call(stack,lapply(1:5479,function(i) setValues(T_min ,round(runif(n = ncell(T_min ),min = -5,max = 24)))))
names(Deu_df_min ) <- paste0('X',gsub('-','.',ymd('1980.01.01') + days(1:5479)))

Here is what I tried:

To do so, I extracted German' shapefile from Eurostat website, here you could see minimal shapefile on the fly (it is perfectly safe, and if you don't want to download, I will provide reproducible shapefile to use). Please take a quick look original minimal shapefile. And I did raster extraction as follow:

library(raster)
shp <- shapefile('eurostat_NUTS3_29-May-18/deu_adm_2006.shp')
deu_extr_mx <- raster::extract(Deu_df_max ,shp)
deu_extr_mn <- raster::extract(Deu_df_min ,shp)
names(deu_extr_mx ) <- shp$NUTS_ID
names(deu_extr_mx ) <- shp$NUTS_NAME
names(deu_extr_mn ) <- shp$NUTS_ID
names(deu_extr_mn ) <- shp$NUTS_NAME

based on the output of deu_extr_mx and deu_extr_mn, I intend to aggregate grid temperature data for each polygon in each year respectively.

Basically, I want to reshape and aggregate T_max and T_min raster grid for each polygon and do merge in plain text ASCII format. I don't have any solid idea how to make the aggregation and reshape data for multi-layers raster grid in R.

Can anyone give me any possible idea or programmatic approach to make this happen in R? How can I reshape and aggregate the raster grid data for the respective polygon and put the merging results in a simple plain text? Any idea to get this done in R more efficiently?

desired output:

I want to aggregate grid temperature observations for each polygon (minimal shapefile on the fly) in each year respectively and put the aggregation result in plain text ASCII format. Is that doable?

Here is the skeleton of my desired aggregation output in plain text ASCII format (I fill it with the random number):

 date   year month day NUTS_ID  NUTS_NAME        lat         long    tmax        tmin
1 1980-1-1 1980    1  1  DE111      Hamm            38.5        -122.5  12.5000     4.1100
2 1980-1-2 1980    1  2  DE112      Ostvorpommern   38.5        -122.5  12.5445     4.2894
3 1980-1-3 1980    1  3  DE113      Aurich          38.5        -122.5  12.5878     4.4574
4 1980-1-4 1980    1  4  DE114      Ludwigslust     38.5        -122.5  12.6298     4.6144
5 1980-1-5 1980    1  5  DE115      Wilhelmshaven   38.5        -122.5  12.6706     4.7604
6 1980-1-6 1980    1  6  DE116      Mecklenburg     38.5        -122.5  12.7103     4.8959

Is that doable to reshape, merge and aggregate temperature grid in multi-layers raster grid for each polygons in each year in R? Is there any efficient way to make this happen? Any idea?


texts: tap on points on touch screen with openlayers

texts:
I have openlayers 2.13 along with geoserver 2.1.3.

I have a vector layer where points created, after a CQL filter is activated.

I have a click event for those points. In pc/laptop works fine, but in mobiles (iphone4s safari and android vimicro firefox) if I tap with my finger a point, nothing happens.

I only load the openlayers library <script type='text/javascript' src='OpenLayers.js'></script> and zoom/pan on the map works in mobiles. CQL filter also works and the points are render in mobiles.

But I can not tap them.

How can I fix this ?

Thanks in advance

the code (I can provide more, if you think is neccessary):

mylayer = new OpenLayers.Layer.Vector("my layer", {renderers: ["Canvas", "SVG", "VML"]})

//Add a select feature control
var select_feature_control = new OpenLayers.Control.SelectFeature(mylayer, {multiple: false,toggle: true});

map.addControl(select_feature_control);

//Activate the control
select_feature_control.activate();

function selected_feature(event){
//do stuff
}

EDIT

I have this now , trying to use Control.GetFeature

//create layer
mylayer = new OpenLayers.Layer.Vector("mylayer", {renderers: ["Canvas", "SVG", "VML"]})

//set protocol
var myprotocol = new OpenLayers.Protocol.WFS({
    url:  "/geoserver/wfs",
    featureType: "mylayer",
    featureNS: "http://www.mysite.gr"
});

//set control
var select_feature_control = new OpenLayers.Control.GetFeature({
  protocol: myprotocol,
  click: true,
  maxFeatures: 1,
  clickTolerance: 1,
  box: false,
  hover: false
});

//add the control
map.addControl(select_feature_control);

//activate the control
select_feature_control.activate();

//register it an connect it with the function
pins.events.register('featureselected', this, selected_feature);

//the function - gets the id by splitting the fid
function selected_feature(event){
var sf = event.feature.fid;
var sfs = sf.split(".");
var sff = sfs[1];
showtext(sff);
}

thats the part of the cql filter

//pass the filter to the layer
var prot =  new OpenLayers.Protocol.WFS({
    url:  "/geoserver/wfs",
    featureType: "mylayer",
    featureNS: "http://www.mysite.gr",
    defaultFilter: filter_final});

var _CallBack = function(resp) {
    pins.addFeatures(resp.features)
}

var response = prot.read({callback: _CallBack});

//refresh to render POIs
mylayer.refresh({force:true});

I can search the layer and see points on the map. But I still can not select in mobile. In pc, I can click a point, but the cursor of the mouse becomes hourglass and nothing happens. Can you spot what is wrong?

Thanks a million


texts: Problem with saving/processing layers in QGIS

texts:
I recently encountered the following problems with a shapefile in QGIS 2.10.1:

when I simply tried to save the layer as a new layer, there is error message: "Feature geometry not imported (OGR error: )".
when I tried to do some processing with the layer, for example, to intersect it with another layer using the "intersection" algorithm from the processing toolbox, it finished without error, but the resulting shapefile does not show up on the map. Actually, in the attribute table there is 0 feature.

The shapefile displays the zip code of a country and it shows correctly on the map. But it seems whatever processing with it leads to an error. Any idea why this happens and how to fix it?

I attached below the link to the shapefile so that you can play with it to see the problem.

shapefile


texts:
We have portal setup at a URL similar to the following through a reverse proxy:

http://<portal url>/arcgis/home

The WebContextURL parameter is set to this:

https://<portal url>/arcgis

Almost everything is working as it should. The only problem we are having now is when we try to create a an App using a template. The app creates fine but the resulting registered url is wrong. The URL ends up being something like this which returns a 404 when you try to go to it

http://<portal url>/app/...

If you change the app URL to this then the app loads:

http://<portal url>/arcgis/app/...

Why isn't it doing this from the get go? Can I change this behavior some where? I haven't found a spot to change this.


texts: Earth engine error: Error: Invalid argument: 'image' must be of type Image

texts:
I'm trying to export an image for some SST temperatures for a given region ( and imported table EEZ) but I'm getting the error;

  Error: Invalid argument: 'image' must be of type Image

Map.addLayer(EEZ, {}, 'EEZ')
Map.addLayer(Chagos, {}, 'Chagos')

//load an image
var sst = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')
.select('sst')
.filterDate('2017-04-01', '2017-04-30')
.filterBounds(EEZ)
.map(function(image){return image.clip(EEZ)}) ;

print('sst', sst)

Map.setCenter(72.25, -6.8, 6); //lat, long, zoom
Map.addLayer (sst, {'min': 23, 'max': 34, 'palette':"0000ff,32cd32,ffff00,ff8c00,ff0000"});

Export.image.toDrive({
image: sst,
description: 'sst Apr 2017',
scale: 500,
region: EEZ
});

After a search I found this post, and changed the export image code to just select the bands I want, as suggested, using the code below. However I still get the same error.

Export.image.toDrive({
image: sst.select('sst'),
description: 'sst Apr 2017',
scale: 500,
region: EEZ
});


texts:
I need to test the performance of the ArcGIS REST API (Version: 10.51) on the map services I have created. The tester within our team is using swagger but it is not implemented in the ArcGIS REST API.

Are there any ESRI tools (external or REST API) to check the performance of the ArcGIS REST API? Where should I look out for them?


texts: Exporting multiple shapefiles into FGDB with ArcPy (Python)?

texts:
Does anyone have any idea why I am getting the error?

I am trying to export multiple shapefiles to a single File Geodatabase.

>>> import arcpy
>>> arcpy.env.workspace = r"W:/gis/condition assessment/quad c/quad c shapefile"
>>> in_features = ["Cond_Asses_Ditch_C.shp", "Cond_Assess_Inlet_C.shp", "Cond_Assess_Manhole_C.shp", "Cond_Assess_Outlet_C.shp"]
>>> out_location = r"X:/data.gdb"
>>> arcpy.FeatureClassToFeatureClass_conversion(in_features, out_location)
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\conversion.py", line 1694, in FeatureClassToFeatureClass
    raise e
RuntimeError: Object: Error in executing tool


texts: In a shapefile, is it possible to have multiple rows linked to a same map feature?

texts:
I want to know whether it is possible to have more than one row of the dbf file that comes with a shapefile map linked to the same feature of such shapefile. The reason for that is simple: to have datasets in the long-format with the shapefiles.

In case, as I suppose, it is not possible, I will still give more details below so someone might suggest good alternatives.

Suppose you have individual shapefile maps of countries' administrative divisions (AD), i.e. one shp file per country. Now, you have information that you want to save per countries' AD for multiple years, given that information is available for different years depending on the country.

The first naïve temptation would be to organize the dbf of the shapefile map in the wide-dataset format, similarly to the following:

The problem with such an approach is that, as I mentioned, each country has INFO on different years. Thus, in case we later want or need to merge countries together in one single world-map, this will become a mess of many variables each covering just a few countries. In a trivial, non-spatial dataset, the solution would be to simply adopt a long-format for the dataset, like the following:

That would be the ideal solution, because then all the dbf of countries' shapefiles would have the exact same variables, with same names, so merging multiple countries would become easy and intuitive. However, for that one would need more than one row of the dbf pointing to the same map feature (in this case, to the same AD of countries).

That is why I ask:

1) is such a thing possible somehow, i.e. to create shapefile maps that have more than one row linked to the same feature of the same map?

If not at all possible:

2) are there any alternatives? I have considered geodatabases but want to avoid file formats that require proprietary format to be read.

My last resource would be to just have two dbfs per country, one with the variables "Country" and "AD" in example I give here, and another with "YR" and "INFO" in long-format, both linked with an ID of some sort. But I am looking for more elegant ways of doing so if any is available.


texts: QGIS Server loads only one .qgs project file

texts: I installed QGIS Server in WinXP via OSGeo4W and runs great. But the problem is, it loads only one .qgs project file. I restarted Apache but it is same. I use QGIS 1.9 default WMS client. I noticed that QGIS Server only load one project file and it looks to the filename first char... For example if i generated  project files named as abc.qgs and xyz.qgs it loads abc.qgs. Do you have any idea ???


texts:
I work in a Geography Department at a small university and we are looking at putting in a 3D Stereo Compilation systems.  We have been looking at Stereo Analyst, PurView and Dat/em, but I was wondering if there are any Open Source alternatives, or a plugin for QGIS.

Any help would be appreciated.


texts:
I'm using MS Visual Studio 2013, C++ and newest gdal downloaded from here http://www.gisinternals.com/query.html?content=filelist&file=release-1800-x64-dev.zip
In readme and libraries I can see HDF 4 and 5, but when I'm doing GDALAllRegister(); and calling

GDALGetDriverShortName(hDriver),
GDALGetDriverLongName(hDriver));

to check drivers, I'm getting only this list of drivers:

As you can see, there is no HDF. How to turn it on or register?


texts: Steps to develop 2.5D web mapping application

texts:
I am developing a 2.5D GIS application. We are using GeoServer for layers.
We created building images which supports 2.5/3D
Already we have developed GIS application using OpenLayers, MapFish etc.

Please any body can guide me to develop this application?


texts: Programatically transfer attributes from line to selected points: : Attribute transfer tool

texts:
Spatial Adjustment>Attribute Transfer Tool in ArcMap transfers attributes between mapped source and target layers by mouse clicks. The disadvantage is you can do one point/line at a time. My problem is to transfer attributes from a single line (a street) to multiple selected points (selected address points).

I was looking into the documentation of IAttribute Transfer:
https://desktop.arcgis.com/en/arcobjects/latest/net/webframe.htm#IAttributeTransfer.htm

I would like to click on the street line and then select multiple address points so that attributes from street will be automatically transferred to selected points- same way as the tool currently works for one point.

 IMxDocument mxDoc = m_app.Document as IMxDocument;
  IMap map = mxDoc.FocusMap;
   if (m_app == null)
    return;
  IEditor editor = m_app.FindExtensionByName("ESRI Object Editor") as IEditor;

  //Cast for Attribute Transfer Type and get IAttributeTransfer.
  IAttributeTransferType attTransType = editor as IAttributeTransferType;
  IAttributeTransfer attTransfer = attTransType.AttributeTransfer;
  //Cast for default settings and return source/target layers.
  IAttributeTransferDefaultSettings attTransDefSettings = attTransfer as
  IAttributeTransferDefaultSettings;
  IFeatureLayer2 sourceFL = ReturnLayer(map, attTransDefSettings.SourceName);
  //If no mapping set, sourceFL will be null (null sourcename) so bail
  if (sourceFL == null)
  {
    System.Windows.Forms.MessageBox.Show
      ("Attribute Transfer field mapping not set");
    return;
  }
  //Get target feature layer.
  IFeatureLayer2 tFLayer = ReturnLayer(map, attTransDefSettings.TargetName);
  //Get the fieldmap for the source layer.
  IFieldMap fieldMap = attTransfer.FindFieldMap(sourceFL.FeatureClass, null);
  //Enumerate through each selected feature (selection from source layer)
  IEnumFeature featEnum = editor.EditSelection as IEnumFeature;
  featEnum.Reset();
  IFeature sourceFeature = featEnum.Next();
  if (sourceFeature == null)
  {
    System.Windows.Forms.MessageBox.Show("Please Selecet a Source Feature");
  }
  editor.StartOperation();
  while (sourceFeature != null)

  {
//TODO: Click on source line and then select multiple target points, which transfers attributes

    try
    {
      if (targetFeature != null)
      {
        attTransfer.Transfer(fieldMap, sourceFeature, targetFeature, out success);
      }
    }
    catch (System.Exception ex)
    {
      MessageBox.Show(ex.Message);
    }
    sourceFeature = featEnum.Next();
  }
  editor.StopOperation("Attribute Transfer");
}

However, that code is looking into automated transfer of overlaying lines. I wanted to interactively select source line and target points by clicking and transferring just as Attribute Transfer tool does (in addition to do multiple selection of points) Any pointers or sample codes to achieve this? I am really stuck.


texts: OpenLayers export map to PNG/PDF example returns file with blank canvas

texts:
I built a basic webmap interface with some vector layers and basemaps (eg. OpenStreetMap) in OpenLayers. I tried the map export examples from OL (sources below), scripts seemingly work perfect and exports files in the correct format, but files contain a blank canvas. I did my research but can't find any appropriate solution to my problem. I already got a suggestion it has something to do with the document.querySelectorAll('ol.layer canvas'), but don't know how to solve it. Below the most relevant part of the code.
Source:

Export map PNG: https://openlayers.org/en/latest/examples/export-map.html
Export map PDF: https://openlayers.org/en/latest/examples/export-pdf.html

const rasterLayer = new ol.layer.Tile({
  source: new ol.source.OSM(),
});

const centerCoords = [4.34878, 50.85045];

const startView = new ol.View({
  center: ol.proj.transform(centerCoords, 'EPSG:4326', 'EPSG:3857'),
  zoom: 8.7
});

const map = new ol.Map({
  layers: [rasterLayer],
  target: document.getElementById('map'),
  view: startView
});

document.getElementById('export-png').addEventListener('click', function () {
  map.once('rendercomplete', function () {
    var mapCanvas = document.createElement('canvas');
    var size = map.getSize();
    mapCanvas.width = size[0];
    mapCanvas.height = size[1];
    var mapContext = mapCanvas.getContext('2d');
    Array.prototype.forEach.call(
      document.querySelectorAll('ol.layer canvas'),
      function (canvas) {
        if (canvas.width > 0) {
          var opacity = canvas.parentNode.style.opacity;
          mapContext.globalAlpha = opacity === '' ? 1 : Number(opacity);
          var transform = canvas.style.transform;
          // Get the transform parameters from the style's transform matrix
          var matrix = transform
            .match(/^matrix\(([^\(]*)\)$/)[1]
            .split(',')
            .map(Number);
          // Apply the transform to the export map context
          CanvasRenderingContext2D.prototype.setTransform.apply(
            mapContext,
            matrix
          );
          mapContext.drawImage(canvas, 0, 0);
        }
      }
    );
    if (navigator.msSaveBlob) {
      // link download attribuute does not work on MS browsers
      navigator.msSaveBlob(mapCanvas.msToBlob(), 'map.png');
    } else {
      var link = document.getElementById('image-download');
      link.href = mapCanvas.toDataURL();
      link.click();
    }
  });
  map.renderSync();
});


texts:
At the risk of a series of down votes but I have just installed QGIS 2.0 after using 1.7.4 and 1.8 previous.

It looks great but there are a lot of plugins missing.  I understand from the install that  I may need to install them again for fear of system conflicts.  However some really good plugins are missing. The following are the main ones:

Numerical Vertex Edit
Easy custom Labelling
OpenLayers Plugin
OpenStreetMap
QGIS Cloud
Remove Empty Layers
Spatial Query

This list is not exhaustive but I am sure they are built in somewhere and can't find them.


texts: How to SELECT all points within a GeoJSON polygon in PostGIS

texts:
I have a GeoJSON polygon, and a table with lots of points. I want to SELECT all points that are within an arbitrary GeoJSON polygon.

How can I do that with PostGIS SQL?


texts: Trying to map student locations and need reference data

texts:
I have a dataset with addresses that I need to map of students. The dataset has the complete address. I am just trying to show pins of where these students exist on a map.
I am trying to create an address locator in ArcGIS (Desktop and Pro) (which is what I think I need but I could be wrong). But to create an address locator, you need a reference file. I don't see anything in help on where to obtain these files.
I tried to download some TIGer files from the census but the ones I found didn't seem to have the required fields I needed for my address search.
Where do you get these reference files?


texts: Seeking shapefile of urban, suburban, and rural areas of US

texts: I'm looking for something like this, doesn't have to be super accurate. I've tried searching for the source of that image but to no avail.


texts: How do you convert a SOAP PolylineN from one namespace to another?

texts:
I'm doing a proof-of-concept to try to speed up a data service that accesses DB2 databases and a GeometryService on an ArcGIS server using the ADF and the SOAP API. I'm trying to speed up data access by getting features using a feature service.

The problem I'm running into is that ArcGIS Server's ADF doesn't seem to be able to convert its own data types.
I have a PolylineN value object from a FeatureService that I need to convert to the same class for a GeometryService.

The documentation online says to use .NET XML serialization, but the examples only convert a service's value object to an ADF value object, not from ADF to a service's value object or one service to another.

Converting from one service's PolylineN to the other doesn't work. I tried a two-step process, converting from FeatureService.PolylineN to ADF.PolyLineN, then from ADF.PolylineN to GeometryService.PolylineN. On deserializing to GeometryService.PolylineN, I get this error:

System.InvalidOperationException: There is an error in XML document (9, 10).
The Inner Exception is the same type, with the message: "The specified type was not recognized: name='PointN', namespace='http://www.esri.com/schemas/ArcGIS/10.0', at ."

According to the documentation, PolylineN is essentially identical in all of these namespaces. So why can't they understand each other's schema?

My workaround is to build a new GeometryService.PolylineN from the feature polyline. It works, but it seems unnecessary to me.

Does anybody have any hints on what I'm missing here?


texts: Read raster at request size through Python

texts:
I'm trying to read the raster (tif) file using Python.
What I want is to read the raster at a specific request size rather than its original resolution.

Example: If a raster has a cellsize of 5 x 5, I would want to read the raster at a lower resolution of suppose 50 x 50 which would reduce the total number of cols and rows. Basically read the pyramid of the raster which has the closet number cols to that of the request size.

Is there a way to read the raster at a particular request size specified through the ArcPy module or gdal module?


texts: Is there any tool for coordinate system detection from a point coordinate?

texts:
I have bounding box points below;

  -20037508.34, -20037508.34, 20037508.34, 20037508.34

Let's think about a tool that gives one of these points as an input and returns the coordinate system/projection which is appropriate. Is there any tool to do this job?

  "-20037508.34" --> Tool -> ED50

By the way what is the projection of coordinates showed above?


texts: How to generate random points that lie on a line in QGIS?

texts: I have a shapefile with track lines and I need to generate random points from each 1 km of every line. There is a function in QGIS Random points but it offers to generate random points only inside a polygon. Any ideas for a specific plug-in or another solution?


texts:
I have to store some polygons in a MySQL database. The polygons are in Plate Carrée (ESPG 32662) projection, but I also tried to first converting them to the usual latitude-longitude projection (ESPG 4326) and it did not work either.

First I upload the shapely polygon to the database:

polygon_text = my_shapely_polygon.to_wkt()
query = f"""
                INSERT INTO
                    lands (polygon_id, geometry)
                VALUES
                    ({id}, ST_GeomFromText({polygon_text}));
        """

The type of the column geometry in the database is set to POLYGON.
However, when I retrieve the geometry and do

shapely.wkb.loads(retrieved_polygon)

I get a point, not a polygon.
The same thing happens if I try to upload the polygon as WKB.


texts: Convert a multiband raster to a polygon in qgis

texts:
My version of QGIS is 2.8.

I want to convert a RGB raster into a polygon. I just need the boundary of it (not the rectangle extent but the boundary of the real feature). I have tried raster calculator but it says SAGA does not support multiband raster.

I also have tried the Vectorize raster layer tool but it takes a very long time because the resolution is very high. It gives me around 3000000 polygon for an image like below and I have a few hundreds this size of images.

Is there any other way to do it?


texts: Creating a bounding box to pull from cache

texts: I have a simple python script that is pulling in lat,longs from company's db and we are trying to create several bounding boxes around these points to simulate zoom-levels. The script will pull tiles from the geoserver, but will not pull from cache. In the geoserver url, I have tiling enabled, but I have been told that the BB's I created do not line up with the "grid" for tiling, so it won't pull from geoservers cach? The documentation makes it seem like i should be able to pull from cache anyway.....any ideas?


texts: Creating map projection for net of geodesic dome?

texts:
I am working on a project that involves mapping the night sky's northern hemisphere onto the exterior surface of a small (16 cm radius) geodesic dome. The problem is, I don't know what map projection to convert to, or whether I need to create one. If so, I don't really know how.

I'm starting from a high-resolution image of the full night sky in an equirectangular projection in celestial coordinates.

I would like to map that onto the net of a geodesic dome,  so that I can construct the dome with accurate star coordinates. I would like to use a computer program to convert the northern hemisphere of my equirectangular star map into the desired projection, if possible.


texts:
im currently building some android application with arcgis, but my client using Identity Manager for auth user on the portal,

i have tried using oauth2 Rest with token in my local development and going well, but when using client portal it's not working,

is there Identity Manager ArcGIS API for Android sample or some kind of documentation on the grid that i can read ?


texts: How to track the the duration of sending data to Google Drive/ Google cloud storage from Google Earth Engine Python API

texts:
I am using Google Earth Engine Python API to claculate NDVI and Export the result in Google Drive / Google Cloud storage (bucket). I am using the following code to export data to Cloud

task = Export.map.toCloudStorage(
    image = toexport,
    #image = SAVI.visualize({min:0, max:255}),
    description = 'Automation',
    scale = 20,
    bucket = 'gisagnext',
    #fileNamePrefix = 'Automation',
    #path = 'savi/',
    path = 'clip/',
    #path = 'ndvi',
    fileFormat = 'png',
    writePublicTiles = 'true',
    maxZoom = 14,
    minZoom = 4,
    #scale = 20,
    #region = row
    region = geometry1
    #maxPixcels=76904329300
    )
task.start()

Now I want to track back the exact time duration of sending data to Google Cloud in Python API.


texts:
I am getting an error when I try to change the value of a feature's Rotation Value.
I have a schematic feature that I would like to rotate by changing the field's ROTATION value, but everytime I try to set the new value I get this error.
"Error HRESULT E_FAIL has been returned from a call to a COM component"

Here is the part of the code where I set the value

ISchematicRelationController controller = new SchematicRelationControllerClass();
var parentContainer = controller.FindParent(duct);
if (parentContainer == null) {
   var rotIdx = parentContainer.Fields.FindField("ROTATION");
   var idx = duct.SchematicElement.Fields.FindField("ROTATION");
   if (rotIdx > -1) {
      Duct.SchematicElement.Value[rotIdx] = 45;
      parentContainer.Value[rotIdx] = 65; ------> this is where the exception occurs
                    }


texts: Determining angle for polyline geometries that share a point

texts:
With help I've written code to determine the angle between two line geometries. However, the code assumes each line extends into infinity and doesn't account for the shared point. Thus an angle of 150 can be returned as 30, or 20 can be returned as 160. I can't come up with a way of refining my code to return the proper angle.

Here's the code:

def AngleFromLines (lines):
    #lines is a python list of line geometries that share a vertex
    for line1 in lines:
        for line2 in lines:
            if line1 == line2:
                continue
            line1StPnt, line1EndPnt = LineToXYs (line1) #get start and end xys for first line
            line2StPnt, line2EndPnt  = LineToXYs (line2)#get start and end xys for second line
            angle1 = GetAngle (line1StPnt, line1EndPnt) #calculate declination
            angle2 = GetAngle (line2StPnt, line2EndPnt) #calculate declination
            angle = abs (angle1 - angle2) #angle between points - doesn't take into account shared point
            #further code not shown....

def LineToXYs (line): #return first and last coordinates
    firstX, firstY = (line.firstPoint.X, line.firstPoint.Y)
    lastX, lastY = (line.lastPoint.X, line.lastPoint.Y)
    return [(firstX, firstY), (lastX, lastY)]

def GetAngle (p1, p2):
    x1, y1 = p1
    x2, y2 = p2
    dX = x2 - x1
    dY = y2 - y1
    rads = math.atan2 (-dY, dX)
    return math.degrees (rads)

AngleFromLines ([lineGeom1, lineGeom2]) #lineGeom* = arcpy polyline geometry object

Example, where the angle is calculated as 39 instead of the desired result of 141.


texts:
I've got a little problem with OpenLayers.
I've a simple function which sets the center of the map to the position of the user - this works fine. But: it should also zoom - which, for an unknown reason, don't work.

Any ideas?

function gotSomeLocation(position){
    var lat = position.coords.latitude;
    var long = position.coords.longitude;
    map.zoomTo(18);
    alert("Lokalisiert!");
    map.setCenter(tCoord(long, lat));

}


texts: Programmatically change layer position in the table of contents (QGIS)

texts:
I am trying move/position a layer to the bottom of the TOC by using addTopLevelItem on all other layers. I unsure how to reference the layers as QTreeWidgetItem's.

def bg_bottom_level_item(self):
    layers = self.canvas.layers()
    layerlist = []
    for layer in layers:
        layerID = str(layer.id())
        if "Auto_Background" not in layerID:
            layerlist.append(layer)
        else:
            pass
    for layer in layerlist:
        leg = qgis.utils.iface.mainWindow().findChild(QTreeWidget, 'theMapLegend')
        item = QTreeWidgetItem(layer)
        leg.addTopLevelItem(item)

Error:

QTreeWidgetItem(QTreeWidgetItem): argument 1 has unexpected type 'QgsRasterLayer'

Or does anyone have another solution for moving layers around in the TOC?


texts: Can markers changes undo through history/revert to original?

texts:
I'm a contributor to a food foraging map and found that many markers are in the ocean. I then also noticed that I had an undo at the top of gmaps and realised I'd moved 5 trees trying to get to my area.

Is it possible to see marker/map history at maps.google.com? Is it possible to revert markers to original position?

I'm completely open to programming this using the api if necessary - surely this information is somewhere...


texts: Using weight in geometric network to Find Path Upstream

texts:
I have a problem about using weights in geometric network in ArcGIS Desktop 10.2. In the example below, I have a simple network with a source & some junctions. I'm trying to use the Find Path Upstream function to find path from point A to the source. I've created a weight based on the length of the edge feature, but it seem doesn't work. Can anyone tell me or show me:

How can I set up a weight to Find Path Upstream from point A to the source through the shortest path (which means A-B-C-D-E-Source in this example)?
Which weight should I use & how can I set it up if I want to find path from A-F-Source?

@Hornbydd: Of course I did set the flow direction & a weight by length but it can't return the result I expect. You can see in the image below is what I receive when I use the Find Path Upstream from A.


texts: How can I create a map of both Californias?

texts:
I want to create a map of "UniCalifornia" which would show the U.S. state of California and Baja California combined into "One California."

The end result should be an image file that does not show bordering areas (no Arizona, etc.) except for perhaps the bodies of waters (Pacific Ocean to the west and Sea of Cortez/Gulf of California to the east of Baja).

What approach need I take to create this? BTW, I am not a cartographer.

UPDATE

If Steven Kay lets me use his image, I came up with a different design than what I describe above:

The combined Californias, if you squint and turn your head sideways, look like a stylized hand with index finger pointing toward something - the Equator, maybe?


texts:
Regarding installing (build) pgrouting

I could install pgrouting-2.0.0 by using a command brew install pgrouting
and following this Build manual: http://docs.pgrouting.org/2.0/en/doc/src/installation/build.html

But unfortunately when I tried to type the make command, following error occurred:

~$make
Scanning dependencies of target apsp_warshall
[  2%] Building C object src/apsp_warshall/src/CMakeFiles/apsp_warshall.dir/apsp.c.o
/Users/username/Build/cpp/pgrouting-git/src/apsp_warshall/src/apsp.c:14:10: fatal error: 'executor/spi.h' file not found
#include "executor/spi.h"
         ^
1 error generated.
make[2]: *** [src/apsp_warshall/src/CMakeFiles/apsp_warshall.dir/apsp.c.o] Error 1
make[1]: *** [src/apsp_warshall/src/CMakeFiles/apsp_warshall.dir/all] Error 2
make: *** [all] Error 2

I couldn't figure out why this is happened and have any solutions.


texts: Is there a way to export and import again all Map at once?

texts:
I have a OpenLayer3 map with some arbitrary layers I can choose at runtime. Layers may come from external WMS services or internal PostgreSQL GIS tables or even from drawn features.

Now, I need a way to export all this stuff to reproduce this same situation again later (on same server). Including zoom level and viewport.

Is there a way to do this (KML, GeoJSON, etc)?


texts: Create a view based on table with GeometryCollection column for use in QGIS

texts:
I have a table with a geometry collection (polygon and point).  QGIS can not read the geometry collection, so I want to create a view for the polygon type and the point type. My SQL function for the polygon type looks like this:

CREATE VIEW parcel_polygons AS (
SELECT
row_number() OVER(ORDER BY inspire_local DESC) AS OID
, inspire_local
, inspire_namespace
, label
, nationalcadref
, areavalue
, validfrom
, beginlifespanversion
, zoning
, geom::geometry(Polygon, 25833)
FROM cp_parcel_postgis
WHERE GeometryType(geom) = 'GEOMETRYCOLLECTION(POLYGON)'

I do not get any errors and the view is created. Unfortunately, the view - table is empty.

What am I doing wrong?

ST_Dump worked. However, I now have the problem that in my geom column of the view geometry types are polygon and point and my ID column (OID) doubles the digits. If I could delete the lines, the problem would be solved. But how do I delete lines in a view?

Unfortunately, the following creation of a rule did not succeed:

CREATE RULE point_delete AS ON DELETE TO parcel_polygons_view DO INSTEAD (
  DELETE FROM parcel_polygons_view WHERE GeometryType (geom) = 'POINT';
);


texts: QGIS: Map theme with layer visibility AND style, does it exist?

texts:
I'm using QGIS to do urban visualisation on large scale territory, I'm used to working with Archicad which has an amazing option called "Graphic Overrides", it's a way to create different layer style (for multiples layers) and quickly switch between them.
It's really great to create different visualisation.
I'm now working with QGIS as I'm working with more data, and I haven't find a good workflow to quickly change styles on multiples layers. I'm talking about unique style per layer saved as a "combination", so I can quickly switch between different visualisation.
Map theme is almost that, but it only switch layer visibility. I know I can save layer style independently and apply theme, but when you want to switch multiple layers, the workflow becomes really cumbersome.
What I wish:

have a list with different combination, which saved styles per layer.
being able to switch between them easily, like the map theme.

What is the closest thing to that in QGIS 3.1x ?


texts: Adding source and target fields in attribute table

texts:
My issue is I want to export a CSV-file from QGIS 3.6.2 to Gephi (network analysis software) the layer I want to export is the edges layer for Gephi. However, when importing into Gephi I am missing the "source" and "target" fields in the CSV file.

More specific I have a layer 'trails' with thousand of features on it and want to see which water hole they go to are they connected? but I don't know what the next step is.


texts: Label polygons in an atlas with similar attributes to the mask layer

texts:
I have two similar vector layers, with the same attributes. One is multipart and I used it as mask for the atlas, the other one is single parts and I want to use it to label the polygons inside each multipart polygon.

I have tried with the rule based labeling with this filter

within($geometry, @atlas_geometry)

But it only works sometimes.

This works

This one not.

I tried to check all the label condictions, and I tried to do it with centroids to ensure that the point was inside the polygon of the mask, but no way to make it work.


texts: Adding DEM Band to existing 3-band TIFF file using QGIS?

texts:
I'm quite new in the GIS business and I do not figure out how to do the following thing:

I have two geoTIFF Files, one is a 3+1 bands file (RGB+alpha), the second one is a DEM geoTIFF file containing one band with elevation data.

For a certain application, what I need is a 1+3+1 bands geoTIFF, containing the elevation Data on the first band, and the optical data plus the alpha info in the 2nd to 4th bands.

How can I join or add bands from different geoTIFF's using QGIS, for instance?

I have found a similar question as Adding band to existing GeoTiff using GDAL?, but there does not seem to exist a definitive answer there.


texts:
I am unsure of the best approach to this problem. I can think of some ways of doing it which I will list below but I am looking for the best and most efficient way of doing this.

Problem

Given a variable number of points, create a "route" or line between the points.
Create a bounding box (of 50km[arbitrary]) around the route.
Return all points within this bounding box
Points must be returned in order along the route (A->D)

Note: The radius does not have to be circles. A,B,C,D can use rectangular bounding boxes.

Data

Table with points (x,y - technically lat/lon)
Input points to create a route
Distance from route to include points

We are using PostgreSQL (8.4) and PostGIS (1.3.6) and Python.

This is the solution I came up with. Thoughts? Ideas? Input?

Possible Solution

Create three polygon "tubes" (A->B, B->C, C->D)
Filter points to only those in the polygons.
Subtract B->C from A->B so the polygons do not overlap (no duplicates)
Subtract C->D from B->C so the polygons do not overlap (no duplicates)
For each polygon
Calculate vector (A->B, B->C, C->D)
Find furthest point opposite of vector in the polygon.
Calculate distance from furthest point to each site in polygon.
Return points in order of distance (smallest -> largest)
Join all the points together (A->B points, B->C points, C->D points)

Solution

My approach was naive, after taking some time and looking at the PostGIS methods I came up with this single SQL call. Note its specific to our database but might help someone in the future. Im sure it could be improved as well.

    SELECT aerodrome_id
    FROM
      geo_aerodromecustom_view,
      (SELECT ST_Buffer(
          ST_Transform(
              (SELECT ST_MakeLine(geom)
              FROM
                  (SELECT geom
                  FROM geo_aerodromecustom_view
                  WHERE aerodrome_id IN ('CYOW','CYVR')
                  ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END
                  ) sites
              ),
              900913
          ),
          50000) as geo
      ) as line
    WHERE ST_Intersects(ST_Transform(geom, 900913), line.geo)
    ORDER BY ST_Line_Locate_Point(
      ST_Transform(
          (SELECT ST_MakeLine(geom)
           FROM
              (SELECT geom
               FROM geo_aerodromecustom_view
               WHERE aerodrome_id IN ('CYOW','CYVR')
               ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END
              ) sites
          ),
          900913
      ),
      ST_Transform(geom, 900913)
    )

I had to project the data because I am using PostGIS 1.3.4 which doesnt support Geography type.

Basically what I am doing is I am  using ST_MakeLine and a query to locate "aerodromes" and return their geometry.

I had to order them (using the CASE directive) so that the line would be connected in the right order.

I then project and buffer this line to create a Polygon that I can then use to see what other aerodromes intersect with the buffered polygon.

Using the unbuffered line (route) and the call ST_Line_Locate_Line I then order the discovered aerodromes as they appear appear along the path.


texts: ArcGIS JavaScript API - Draw without Using API Toolbar Widget

texts:
ArcGIS API for JavaScript is using ESRI Toolbar Widget to handle Drawing like Point on the maps. Here is the Demo doing this job.

Now my question is, is there a way to not use the "esri/toolbars/draw" facilities (Widget) and draw on the Graphic Layer by using core API coding?


texts: Cut points falling in a polygon out of the point cloud

texts:
Using lidR 3.1.3, is it possible to classify points inside polygons shapefile for a LASCatalog, something like the previous lasclassify() function. I'm looking for a equivalent of this PDAL method.
ctg <- lidR::readLAScatalog(path_to_folder)
shapefile = rgdal::readOGR(shapefile_file)
opt_output_files(ctg) <- paste0(path_to_folder, "/{ORIGINALFILENAME}_classified")

# Need to classify here
output<- classify_within_polygons(ctg, shapefile)


texts: pgpointcloud repackaging point or x, y, z geometries to patch

texts:
I have a pointcloud in PostgreSQL/postgis imported via PDAL, everything is fine, I chipped it into groups of 600 points, no problems.

I then manipulated the data, first exploding the patches and saving the coordinates to a new table:

CREATE TABLE pointcloud52 AS
SELECT
st_x(PC_EXPLODE(pa)::geometry) as x,
st_y(PC_EXPLODE(pa)::geometry) as y,
st_z(PC_EXPLODE(pa)::geometry) as z
From pointclouds

Now I have columns x,y,z how do I re-patch these, i.e. return the points to a patched geometry?

I've been trying with the PC_Patch() function, but it refuses to work


texts: Is there a method to convert the text of labels in lines? I think infact that the letters can be hypothetically transformed in lines, for simplify the view in other kinds of gis.


texts: Different styles of area in same shape layer

texts:
I'm using QGIS 3.16.
I have the following problems: One of my shape-layer has several areas. Now I want to give every area a different style/colour. I already found the function to give all areas the same style/colour. I had an idea today in the morning that I could multiply or copy the shape-layer and edit in every layer one area and delete the rest. That would make a lot of effort. Does somebody has an other idea to change the areas in the same layer?
In addition, I don't find the function which you can change the graduation of the scale bar. So no 100, 200, 300, 400 steps, instead 50 100 200 400 etc.


texts: Experience with Geoserver's GML Extension

texts: I'd like to know if Geoserver's GML extension works (well) with GeoServer 2.1-RC2 and whether it will enable uploading of GML files via REST API.


texts: How to link the residential buildings and residential in osm?

texts:
In osm. I can map building  with tag building = residential
also can map Residential way with tag landuse = residential
but how to map The building belongs to residential?
I thought of a way  map the building with operator = 'xxxx'. But it does not seem right


texts: Is there a plugins in QGIS that can auto rectify Openstreetmap vector area change sets as per aerial imagery


texts: Parallel Offset of a Polygon Edge with Shapely

texts:
I need to offset edges of a polygon by different amounts.

buffer() will parallel-offset all the edges by the same amount, and each in its edge-normal direction (and in both directions if applied to a LineString).

parallel_offset() will offset a whole LineString by the same amount, but that means that if given only some edges of a polygon, the result is disconnected from the other neighboring polygon edges which were not offset. It will take care of the corners/joints only within the LineString and not on the sides.

Is there a way to nicely offset polygon edges with different amounts while still getting nice correct corners?


texts:
I've been going through the GeoTools tutorials and have read the Vector Grids page http://docs.geotools.org/latest/userguide/extension/grid.html

I am trying to replicate the maps shown on this page with the following code, however I'm unable to generate any output (the program doesn't open the map window or output any errors). Where am I going wrong?

import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.grid.Envelopes;
import org.geotools.grid.GridFeatureBuilder;
import org.geotools.grid.Grids;
import org.geotools.helpers.IntersectionBuilder;
import org.geotools.map.FeatureLayer;
import org.geotools.map.Layer;
import org.geotools.map.MapContent;
import org.geotools.styling.SLD;
import org.geotools.styling.Style;
import org.geotools.swing.JMapFrame;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.geometry.coordinate.Polygon;

import java.io.File;
import java.net.URL;

public class Grid1 {

  public static void main(String[] args) throws Exception {

    // Load the outline of Australia from a shapefile
    // display a data store file chooser dialog for shapefiles
    File file = JFileDataStoreChooser.showOpenFile("shp", null);
    if (file == null) {
      return;
    }

    FileDataStore dataStore = FileDataStoreFinder.getDataStore(file);
    SimpleFeatureSource ozMapSource = dataStore.getFeatureSource();

    // Set the grid size (1 degree) and create a bounding envelope
    // that is neatly aligned with the grid size
    double sideLen = 1.0;
    ReferencedEnvelope gridBounds = Envelopes.expandToInclude(ozMapSource.getBounds(), sideLen);

    // Create a feature type
    SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
    tb.setName("grid");
    tb.add(
        GridFeatureBuilder.DEFAULT_GEOMETRY_ATTRIBUTE_NAME,
        Polygon.class,
        gridBounds.getCoordinateReferenceSystem());
    tb.add("id", Integer.class);
    SimpleFeatureType TYPE = tb.buildFeatureType();

    // Build the grid the custom feature builder class
    GridFeatureBuilder builder = new IntersectionBuilder(TYPE, ozMapSource);
    SimpleFeatureSource grid = Grids.createHexagonalGrid(gridBounds, sideLen, -1, builder);

    // Create a map content and add our shapefile to it
    MapContent map = new MapContent();
    map.setTitle("Quickstart");

    Style style = SLD.createSimpleStyle(ozMapSource.getSchema());
    Layer layer = new FeatureLayer(ozMapSource, style);
    map.addLayer(layer);

    Style layerStyle = SLD.createSimpleStyle(grid.getSchema());
    Layer layerGrid = new FeatureLayer(grid, layerStyle);
    map.addLayer(layerGrid);

    // Now display the map
    JMapFrame.showMap(map);

  }
}


texts: Is there a publicly available dataset that I can use that links 9 digit zipcodes with counties and states?

texts: I have a dataset which has 9 digit zipcodes and I am trying to match these to the state and county FIPS codes. Is there a publicly available dataset that will allow me to do this?


texts: Hydrology tool dataset in ArcGIS 10.1 for calculating slope?

texts:
I am quite new to the hydrology toolset in ArcGIS and thought it may be useful for extracting slope data from a raster. I was wondering if anyone had experience trying to extract slope data upslope from a point?

I basically need to find a location 10 meters (on the ground, not as the crow flies) up the steepest slope from a known X,Y location. As a 10m buffer will not take into account surface geomorphology, I was wondering if anyone had any ideas regarding tools that could be useful?  Seems like drainage tools could be handy in this case, but I'm not sure which ones to use or how to create a polyline from drainage raster datasets.

This would eventually be put into a python script to automate slope extraction for landslides.


texts: Problem with mouse zoom in Open Layers 6

texts:
I'm using latest version of OpenLayers 6 and I noticed that even in the official examples there is small problem (or it is a feature) with mouse zoom. When I set zoom restrictions then while using the zoom in/zoom out buttons I can't zoom more or zoom out according those restrictions. But with mouse wheel it behaves strange, the map is zoomed out (despite restrictions) and when you stop it comes back to its place.

Take a look here: https://www.screencast.com/t/P6skAIhl433V

Is there a way to disable it?


texts: Quickly find all polygons that overlap with two or more polygons in Shapely

texts:
I have a set of Polygons (in Shapely) that I wish to find all intersections that are comprised of overlaps of two or more polygons. I have been working based off the code found in this answer. I have a code so far that works in the if 1: statement below, but it takes a very long time to run. The code that is in the if 0: statement, runs quickly, but produces an output with some errors. See below for an example of the error.

from shapely.geometry import Point, LineString, shape, mapping
from shapely.ops import cascaded_union
import fiona
import shapely
from shapely.ops import unary_union, polygonize
from osgeo import ogr
from shapely.strtree import STRtree

with fiona.open("Untitled.shp") as layer:
   subset=layer
   spoof=[]
   for n,pol in enumerate(subset):
      if type(shape(pol['geometry'])) == shapely.geometry.polygon.Polygon:
         spoof.append(shape(pol['geometry']))
      elif type(shape(pol['geometry'])) == shapely.geometry.multipolygon.MultiPolygon:
         print('theres one here')
         tmp=list(shape(pol['geometry']))
         for x in tmp:
            spoof.append(x)

   rings = [LineString(list(pol.exterior.coords)) for pol in spoof]
   union = unary_union(rings)

   tree = STRtree(spoof)
   result = [geom for geom in polygonize(union)]

   final=[]
   if 1:
      for k,item1 in enumerate(result):
         count=0
         for n,item2 in enumerate(spoof):
            if item1.intersection(item2).area >1e-8:
               count=count+1
         if count>=2:
            final.append(item1)

   if 0:
      for k,item in enumerate(result):
         touch=tree.query(item.centroid)
         if len(touch)>=2:
            final.append(item)

Below, the orange represents all the shapes that I would like to check if they are produced (by unary union in Shapely) by two or more overlaps, the cyan represents the output of the if-statement that works (if 1 in above code), and the pink represents the incorrect selected polygons by the if-statement that doesn't work (if 0 in above code). The top left corners of the text boxes in the pictures represent that point being selected, and as you can see the second picture, the pink area only has one polygon from the STEREO_OBSERVATIONS (which is the set of polygons that are being checked for overlap) so it should not be selected. By contrast, the left adjacent cyan square to the pink square highlighted in the third picture has two polygons from STEREO_OBSERVATIONS and as such is correctly selected by the cyan result. The layering of the polygon layers is shown in the first picture.

What is the best way to query the STRtree structure so as to check to intersection or overlap? I have tried to use the centroid in the code, but that hasn't seemed to produce the correct output.

Fig. 1 (Legend/layering)

Fig 2. (Incorrect selection by fast algorithm [if 0 in code] in pink)

Fig 3. (Correct selection by slow algorithm [if 1 in code]


texts: Create lines to connect vertices in QGIS

texts: I have a file that marks street edges and want to connect "crossing" points automatically. Image 1 shows the street edge, Image 2 shows the desired result. Is there a plugin or method I can use to accomplish this without having to scroll through and create a new line at each intersection?


texts: Merging different classes/clusters in one single raster in ArcGIS Desktop?

texts:
I've created a signature file with 8 clusters and performed an unsupervised classification in ArcGIS 9.3. Later I have grouped values (spatial) classes together to form 4 classes I needed (in symbology) but that just "colors" them the same, they're still treated as separate categories. I need to merge them to have only 4 groups/categories.

What to do?


texts: Mapserver: "Premature end of script headers: mapserv,"

texts:
Some of requests to the Mapserver return "Internal server error". When I checked the Apache log I saw that there are several error messages:

GDAL: In GDALDestroy - unloading GDAL shared library., referer: http://...
GDAL: In GDALDestroy - unloading GDAL shared library., referer: http://...
Premature end of script headers: mapserv, referer: http://...

I'm not familiar with the GDAL and I'm not sure how to fix the problem.

Update

Ubuntu 13.04
GDAL 1.9.0, released 2011/12/29
Mapaserver 6.2.1
POSTGIS="1.5.3"
GEOS="3.3.8-CAPI-1.7.8"

The "Premature headers end" resolves in not loading some part of the map, but not the entire map. The strangest think is that when I zoom in some of the missing parts are getting shown.

Update 2:

Here is another error I found:

msProjectRect(): some points failed to reproject, doing internal sampling.

Update 3:

The main problem was in the extent set for the layer in the map file. It was not cover all the polygons which I want to render.


texts: How to filter Search widget suggestions/results to a single US State?

texts:
I have an esri search widget I'm building like so:

const searchWidget = new Search({
    view: view,
    id: "search_widget",
    locationEnabled: true,
    includeDefaultSources: false,
    sources: [
        {
            name: "CA Locator",
            placeholder: "Enter a CA address...",
            countryCode: "US",
            locator: new Locator({
                url: "https://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer"
            }),
            singleLineFieldName: "SingleLine",
            outFields: ["State", "Addr_type"],
            suffix: ", California, USA",
            filter: {where: "State = 'California'"}
        }
    ]
});

I've tried to add a suffix in attempts to force the searchTerm to resolve to California without luck.

How can I use filter with World GeocodeServer?


texts: How to perform a "vertical integration" using FME?

texts:
I am not a GIS person, but I need to understand an aspect of the project I'm working on.

We are taking an ArcSDE layer, and importing it into FME.  I now need to do some checking of the data.  In ArcSDE I did topological checks to find if the layers were "vertically coherent".  How would I do this in FME?

Vertical coherence: to me vertical coherence is that the boundaries of one polygon is within a certain tolerance level of another in a different layer that "overlaps".  The boundary lines are the same.


texts: Is the reflectance required to get the NDVI, for Landsat 8 images?

texts:
I'm a beginner in satellite optical imagery.
I have a landsat 8 image (11 bands downloaded from http://earthexplorer.usgs.gov/) for which I need to calculate the NDVI (Normalized Difference Vegetation Index) and possibly the NDWI (normalized difference water index) as well.

Is it required to get the reflectance from the DN (calibration) before calculating the NDVI?
If so, what's the open-source software usually used for the purpose of this calibration?


texts: Creating KMZ/KML files with NetworkLink and Regions in Java?

texts:
I have already made KML/KMZ files using org.geotools in Java (see below). But now we have got so large files that I need to cut them into regions and use networkLinks. This is mentioned in https://developers.google.com/kml/documentation/kmlreference#networklink .

How do I do that in Java using org.geotools?
(Nothing goes wrong with the code below, but it does not include regions and networkLinks.)

public void export(Collection<Occurrence> data, OutputStream os)
{
  SimpleFeatureCollection features = createSampleFeatures(data, true);
  Encoder encoder = new Encoder(new KMLConfiguration());
  encoder.setIndenting(true);

  try (ByteArrayOutputStream baos = new ByteArrayOutputStream();
    ZipOutputStream zos = new ZipOutputStream(os);
    ByteArrayOutputStream tempOs = new ByteArrayOutputStream())
    {
       encoder.encode(features, KML.kml, tempOs);

       deleteRootElement(baos, tempOs);

       ZipEntry ze = new ZipEntry("doc.kml");

       zos.putNextEntry(ze);

       PrintStream ps = new PrintStream(zos);
       ps.println(baos.toString());
       ps.flush();
       zos.closeEntry();
       }
       catch (IOException | TransformerFactoryConfigurationError | TransformerException e)
       {

       }
    }
}

with

protected SimpleFeatureCollection createSampleFeatures(
   Collection<Occurrence> occurrences, boolean withLink)
{
  SimpleFeatureTypeBuilder typeBuilder = new SimpleFeatureTypeBuilder();
  typeBuilder.setName("occurrences");
  typeBuilder.add("geometry", Geometry.class, DefaultGeographicCRS.WGS84);
  typeBuilder.add("name", String.class);
  typeBuilder.add("description", Object.class);
  SimpleFeatureType type = typeBuilder.buildFeatureType();

  SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(type);
  List<SimpleFeature> simpleFeatures = new ArrayList<>();
  CoordinateTransformerFactory factory = new CoordinateTransformerFactory();

  Iterator<Occurrence> iterator = occurrences.iterator();
  while (iterator.hasNext())
  {
    Occurrence occurrence = iterator.next();
    String wkt = 4326;

    WKTReader wktReader = new WKTReader();
    try
    {
      Geometry geom = wktReader.read(wkt);
      featureBuilder.add(geom);
    }
    catch (ParseException e)
    {
    }
    featureBuilder.add(occurrence.getTaxon().getName());
    String link = "this is a link";

    featureBuilder
       .add(link + " "
         + StringFormatter.getNullSafeStringRepresentation(
           occurrence.getSample().getSurvey().getAvailability())
         + " " + occurrence.getSample().getDate() + " Fund-UUID: "
         + occurrence.getUuid());
       simpleFeatures.add(featureBuilder
       .buildFeature("uuid." + occurrence.getUuid().toString()));
     }

     return DataUtilities.collection(simpleFeatures);
}

and

// Workaround for https://trac.osgeo.org/gdal/ticket/6981
// remove after fix is published for recent gdal and qgis versions

protected void deleteRootElement(OutputStream os, ByteArrayOutputStream tempOs)
 throws TransformerFactoryConfigurationError, TransformerException
 {
   TransformerFactory factory = TransformerFactory.newInstance();
   Source xslt = new StreamSource(this.getClass().getResourceAsStream("/removeNs.xslt"));
  Transformer transformer = factory.newTransformer(xslt);
  Source text = new StreamSource(new ByteArrayInputStream(tempOs.toByteArray()));
  transformer.transform(text, new StreamResult(os));
}


texts: I used the method extract by mask  to clip a raster , but  the output raster beyond the boundary, do you know why ?


texts:
How can I crop an area from an aerial image knowing the latitude and longitude of the image's corners and also the cropped section?

I have a LiDAR scan and I and want to crop an image overlay to fit over the scanned area.


texts: Minimum bounding geometry with "group" option gives error in ArcMap

texts:
The minimum bounding geometry tool is working fine on my data, but when I try using the "group" option, I get ERROR 999999: Error executing function.

I would like generate two polygons, A and B, which are the convex hull of their corresponding group of points. This is what the attribute table of my point shapefile looks like:

lat  long  ID
Ax1  Ay1   nameA
Ax2  Ay2   nameA
Ax3  Ay3   nameA
Bx1  By1   nameB
Bx2  By2   nameB
Bx3  By3   nameB

I can run the MBG (Minimum Bounding Geometry) tool for all the data and get one big polygon. But this isn't what I want; I want individual polygons based on their ID. When I use the "group" function within the MBG tool, I select LIST and my ID column. But the process fails and I get the aforementioned error.

Other details:
I am using ArcMap 10.3.1.
The exact error log is as follows:

Input features

sheet$ Events
Geometry Type: Convex Hull
Group Option: List
Group Field: ID
Add geometry characteristics as attributes to output: NO_MBG_FIELDS

Messages:

  ERROR 999999: Error executing function.

  Failed to execute (MinimumBoundingGeometry)


texts: Exporting lists of Feature classes and fields from SDE

texts:
I'm using matt wilkies' arcplus module https://gis.stackexchange.com/a/5943/16793 to export a list of feature classes within an SDE database to a text file.   I've also written a tool to write field names to a csv file:

   """The source code of the tool."""
  import arcpy, os, csv

fclass = "SDEConnection\\"+Feature
fieldnames = [f.name for f in arcpy.ListFields(fclass)]
writepath = r"C:\Users\ME\Desktop\TableFields\\"+Feature+".csv"
with open(writepath, "w") as output:
        writer = csv.writer(output, lineterminator='\n')
        for val in fieldnames:
    writer.writerow([val])

Ideally, I'd like to loop through the feature classes and fields and write a write file containing something like below:

Feature Class or Table-----------------------------------Fields

WELLS.WELL_DATA-------------------------------------WELL_KEY, SITEID, DEPTH, DIAMETER

I'm a little stumped on nesting the loops and formatting the output like this.

-----------EDIT-----------------------------------------------------------

Using the Walk method with a script like this

import os, arcpy, csv
def inventory_data(workspace, datatypes):
for path, path_names, data_names in arcpy.da.Walk(
        workspace, datatype=datatypes):
    for data_name in data_names:
        yield os.path.join(path, data_name)

obtains a unicode string like "C:\Users\user\AppData\Roaming\Esri\Desktop10.2\ArcCatalog\Connection to oracle.sde\BASINS\BASINS.ONE"

while

fn = [f.name for f in arcpy.ListFields(fc)]
print fn

extracts a list

[u'FNODE_', u'TNODE_', u'LPOLY_', u'RPOLY_',]

My aim is to combine Walk with ListFields in a loop and write the FeatureClass with its accompanying fields (delimited) to a file.