texts: Extending Polygon to Another Polygon with Automation?

texts: I need to extend the polygons in blue to match the red polygons using automation, since there are thousands of these. It seems like a simple thing to do, but I have not found a suitable way to accomplish it.


texts: Wetlands Symbol for QGIS point pattern fill

texts: Is there a marker available anywhere for point pattern fill that resembles the standard wetlands symbol used in USGS topographic and other maps? Any good standard symbol that will be intuitively recognized as wetlands?


texts:
I am trying to insert Field Values to an empty table.

My code by so far:

arcpy.CreateTable_management(r"C:\Users\esri\Documents\ArcGIS","test")
arcpy.AddField_management("test","total","LONG")
arcpy.AddField_management("test","incomplete","LONG")
arcpy.AddField_management("test","completed","LONG")
arcpy.AddField_management("test","notstarted","LONG")
list=[1,2,3,4]
ins=arcpy.da.InsertCursor("test","*")
ins.insertRow(list)

It shows the error, the sequence must be matched. I am also tried with

ins.insertRow([list])


texts: Setting up PyCharm 64bit with QGIS 64bit OSGeo4W64 installation on Windows

texts:
Trying to setup PyCharm Community Edition 64Bit with QGIS OSGeo4W64 installation on Windows so that I can debug some python and reference the QGIS/Qt modules.

I tried the folowing tutorials;
https://planet.qgis.org/planet/tag/pycharm/
Adding QGIS-specific Python interpreter to PyCharm on Windows?
https://silverspringenergy.com/using-pycharm-as-an-ide-for-qgis-3-plugin-development-2/

The last one creates a .bat file that is suppose to automatically open PyCharm but when it runs nothing opens, no errors either.

@echo off
call "C:\OSGEO4W64"\bin\o4w_env.bat
call "C:\OSGEO4W64"\apps\grass\grass78\etc\env.bat

@echo off
path %PATH%;C:\OSGEO4W64\apps\qgis\bin
path %PATH%;C:\OSGEO4W64\apps\grass\grass78\lib
path %PATH%;C:\OSGeo4W64\bin\apps\Qt5\bin
path %PATH%;C:\OSGeo4W64\bin\apps\Python37\Scripts
set PYTHONPATH=%PYTHONPATH%;C:\OSGeo4W64\bin\apps\qgis\python
set PYTHONHOME=C:\OSGeo4W64\bin\apps\Python37

start "PyCharm aware of Quantum GIS" /B "C:\Program Files\JetBrains\PyCharm
Community Edition 2019.3.1\bin\pycharm64.exe"

Could someone please help me setup the Community Edition of PyCharm with the OSGeo QGIS installation on a Windows OS?


texts:
Is there a way to find the compactness (ratio of perimeter to area) of various groups of contiguous cells using GRASS functions?

If you imagine I want to find the compactness of the three land parcels shown in the below picture

Perhaps there are some reporting functions I can use to find the area and perimeter of these parcels and then I can determine the ratio myself?


texts:
I'm using ArcMap 10.1. I want to select all vertices that contain z values and decrease their value with -0,155 meter but I dont want to change vertices that have the value=0.
With this script I managed to change all vertices but how can I change it to only update those with a value?

fc = r'C:\Users\tomlin01\Documents\ArcGIS\Default.gdb\Export_Output_8'    #layer to update
z_increase = -0,155         #change in meters
with arcpy.da.UpdateCursor(fc, 'SHAPE@Z', explode_to_points=True) as cursor:
    for row in cursor:
        cursor.updateRow([row[0] + z_increase])


texts: Using Google Maps API v2 and Maps engine?

texts:
I am researching on the options in drawing shapes on Google Maps with API v2 and Google Maps Engine. I have managed to set up the admission to draw a viewed layer on my setup on my site.

Does anyone know if it is possible to access your own layers and maps and draw an save the tables with attributes and so on?


texts: Addressing fields in ArcGIS field calculator dynamically?

texts:
I am (re-)learning ArcGIS for a new job after couple of years on PostGIS (and QGIS/pyqgis) and can't figure out how I can address a Field within the pre-logic-block of the field calculator.

Pre-Logic-Field:

def createStyleField(kindofuse):
    if kindofuse == 'wa':
        listofrelevantfields= {
            1: ['level1', 'level1_wa'],
            2: ['level2', 'level2_wa']
        }
    elif kindofuse == 'ge':
        listofrelevantfields= {
            1: ['level1', 'level1_ge'],
            2: ['level2', 'level2_ge']
        }
    elif kindofuse == 'gi':
        listofrelevantfields= {
            1: ['level1', 'level1_gi'],
            2: ['level2', 'level2_gi']
        }
    results = {}
    for levelkuerzel, liste in listofrelevantfields.iteritems():
        for entry in liste:
            if '!' + entry + '!' == 1: ## THIS IS WHERE I AM TRYING TO ADDRESS A FIELD
                results[levelkuerzel] = True
                break
            else:
                continue
    for i in range(1, 5) : #up to 5 levels, the lowest will be returned
        try:
            if results[i] == True:
                return unicode(i)
        except KeyError as err:
            continue
        except Exception as err:
            raise err
    return 0

And the other Field:

createStyleField('wa')

I need this for a couple of more complex calculations with a big number and case-dependent selection of fields. How do I do this?

Or should I address this issue totally different?


texts: Decreasing distance between legend's items in R?

texts:
I'm trying to create correct choropleth() map in R using GISTools package. I use choro.legend() to display my legend. However I can't create either correct legend placement, AND/OR I don't understand the role of fmt parameter in my choro.legend() function. It seems to me that fmt can decrease the space between my legend colors and description.

I found choro.legend help page that states fmtis "C style format for values stated in above choropleth class limits".

So this should just work with my classes values, not by the distance between legend items itself?

OR what is the human understandable meaning of this?

How can I decrease the distance between my choro.legend()'s items?

choropleth(my.shp, nc.lI[,1], shading = income.shade)
choro.legend(-12919698, 5314317,income.shade,title='My title',cex=0.8, bty = "n", fmt = "%0.1f")


texts: Algorithms for optimally locating points

texts:
I am trying to compare locations of where several thousand facilities have been actually built to where they would be optimally located to minimize travel times of the population (represented by census block or tract centroids).  I'm having trouble finding much of anything on how to optimally locate points.

I've got some idea of how to go about choosing these locations, but the sheer number of points to be placed in space means that any non-cleverly optimized algorithm is going to take a long time, possibly years.  Thus my question: Are there standard algorithms for choosing where to locate a fixed number of points?

I will ultimately take whatever algorithm I find as a starting point and adapt it to incorporate more information than just the population counts.  Thus the preferred answer would include a detailed description of the algorithm, code, or be written in an open-source language, so that I can replicate and extend it.  However, if ArcGIS has a convenient function for this optimization, I'd be happy to start with that.


texts: Custom Printing Service with Data Driven Pages

texts:
I've been tasked with creating a custom web app hosted in ArcGIS Online that has the ability to export out data driven pages through a printing service (we also have ArcGIS Server).  Within the printing widget, I noticed you can point to a printing service.

What we want to be able to do is have any user query features, go through the printing widget and have it export out a multi-page PDF of all features that are queried.  I am familiar with sharing the "Get Layout Templates Info" as a geoprocessing service but I'm running into a bit of trouble with the data driven pages aspect.


texts: Execute multiple process in parallel of Arcpy script

texts:
I have seen some difficult ways to do parallel processing, but I wonder if it is possible to simply execute multiple process of the same ArcPy script at the same time.

My script makes some changes to the default geodatabase, so I thought of making a geodatabase copy for each process.

I have updated the script to copy the shared resources between the processes, so it copies the geodatabase and the mxd's related to it.

I have made a test to parallelize, using this script:

pool = multiprocessing.Pool(2)

pool.map(test_func, [1, 2] , 1)

pool.close()

I noticed when I browse RAM and CPU use, every process consumes 200Mb. So, if I have 6 Gb of RAM, I think I have to exploit 5Gb RAM of it by enlarging the pool size to:

5000 / 200 = 25

So, to exploit the whole power of the machine, I think I should use 25 as pool size.

I need to know if this is the best manner, or how I could measure the efficiency of this parallelization.

This an example of the code that I'm trying to parallelize. The whole script contains 1500 lines of code almost like this one:

def dora_layer_goned():
    arcpy.Select_analysis( "layer_goned" , "layer_goned22" )
    arcpy.MakeFeatureLayer_management("layer_goned", "layer_goned_lyr")
    arcpy.SelectLayerByLocation_management("layer_goned_lyr" ,"WITHIN", "current_parcel"  ,  "" , "NEW_SELECTION")
    arcpy.SelectLayerByAttribute_management("layer_goned_lyr" , "SWITCH_SELECTION" )
    arcpy.Select_analysis("layer_goned_lyr" , "layer_goned_2_dora2" )
    arcpy.Clip_analysis("layer_goned_lyr" , "current_parcel_5m_2" , "layer_goned_2_dora" )
    arcpy.SelectLayerByAttribute_management("layer_goned_lyr" , "CLEAR_SELECTION" )
    arcpy.FeatureToPoint_management("layer_goned_2_dora","layer_goned_2_dora_point","CENTROID")
    arcpy.MakeFeatureLayer_management("layer_goned_2_dora2", "layer_goned_2_dora_lyr")
    arcpy.SelectLayerByLocation_management("layer_goned_2_dora_lyr" ,"INTERSECT", "layer_goned_2_dora_point"  ,  "" , "NEW_SELECTION")
    arcpy.DeleteFeatures_management("layer_goned_2_dora_lyr")
    arcpy.FeatureVerticesToPoints_management("current_parcel","current_parcel__point", "ALL")
    arcpy.FeatureVerticesToPoints_management("carre_line","carre_line__point", "ALL")
    arcpy.CalculateField_management("current_parcel__point","id","!objectid!","PYTHON_9.3")
    arcpy.SpatialJoin_analysis("carre_line__point"  , "current_parcel__point" , "carre_line__point_sj","JOIN_ONE_TO_ONE" , "KEEP_COMMON" , "" , "CLOSEST")
    arcpy.Append_management("current_parcel__point" , "carre_line__point_sj" , "NO_TEST") #
    arcpy.PointsToLine_management("carre_line__point_sj", "carre_line__point_sj_line", "id")
    arcpy.Buffer_analysis("carre_line__point_sj_line" , "carre_line__point_sj_line_buf" , 0.2)
    arcpy.Erase_analysis("layer_goned_2_dora2" , "carre_line__point_sj_line_buf"  , "layer_goned_2_dora_erz")
    arcpy.MultipartToSinglepart_management("layer_goned_2_dora_erz" , "layer_goned_2_dora_erz_mono")
    arcpy.MakeFeatureLayer_management("layer_goned_2_dora_erz_mono", "layer_goned_2_dora_erz_lyr")
    arcpy.SelectLayerByLocation_management("layer_goned_2_dora_erz_lyr" ,"SHARE_A_LINE_SEGMENT_WITH", "current_parcel"  ,  "" , "NEW_SELECTION")
    arcpy.SelectLayerByLocation_management("layer_goned_lyr" ,"CONTAINS", "layer_goned_2_dora_erz_lyr"  ,  "" , "NEW_SELECTION")
    arcpy.DeleteFeatures_management("layer_goned_lyr")
    arcpy.Append_management("layer_goned_2_dora_erz_lyr" , "layer_goned" , "NO_TEST") #
    arcpy.SelectLayerByAttribute_management("layer_goned_2_dora_erz_lyr" , "CLEAR_SELECTION" )


texts: Attribute from point to line using end Vertex of a line - py console

texts: I have a line feature and a point feature. I would like to join the attribute of a point feature to the line feature based on end vertex, so when end vertex of a line falls on a point, the line feature gets the attribute of the 'point'(As show in the  Pic). Would like to know how to use py console. or is there a script already in place?


texts: X and Y after Albers lat and lon transformation

texts:
that is probably a very obvious question but it is melting my brain:

when transforming lat and lon to [x,y] (for instance using this code https://gist.github.com/RandomEtc/476238) I get the following:
{ x: -0.27, y: -1.12 }

What are x and y? Are these radians? How can I use x and y to find the exact point for that lat and lon in may 2D projected map?

Edit:

This is the map I am working with: http://www.jetplan.com/weather/data/maps/eutball06.gif

What I am trying to achieve is to map a (lat,lon) to a (x,y) pixel on that image.

I think the map is an Albers projection.

I have used the code mentioned above to get an x and y value.

My problem is that the (x,y) coordinates I get are related to a unit sphere and I am not sure how am I supposed to map this to a pixel on the image of the map.


texts: `ogr2ogr` command doesn't run the same through Docker

texts:
I have the following command:

ogr2ogr --debug ON -f 'KML' ./areas_desarrollo_indigena_2017_03_07.kml PG:"host=localhost port=5432 dbname=mapserver user=postgres password=" -sql "SELECT * from ( SELECT ST_Intersection(datasrc.geom, polygon) as sharedPolygon, datasrc.* from ST_MakePolygon(ST_GeomFromText('LINESTRING(-69.02160644531251 -20.053351219365318,-69.78515625000001 -20.45789615504684,-69.09851074218751 -20.979392255760608,-69.02160644531251 -20.053351219365318)',4326)) as polygon,datasources.areas_desarrollo_indigena_2017_03_07_15 as datasrc) as sharedPolygon WHERE ST_IsEmpty(sharedPolygon)='false'"

Which runs perfectly. Notice the added --debug ON, I've run it with and without debug.

Now I want to execute the same command through docker because I want to target a specific GDAL version:

docker run "geodata/gdal:2.1.3" ogr2ogr --debug ON -f 'KML' ./areas_desarrollo_indigena_2017_03_07.kml PG:"host=localhost port=5432 dbname=mapserver user=postgres password=" -sql "SELECT * from ( SELECT ST_Intersection(datasrc.geom, polygon) as sharedPolygon, datasrc.* from ST_MakePolygon(ST_GeomFromText('LINESTRING(-69.02160644531251 -20.053351219365318,-69.78515625000001 -20.45789615504684,-69.09851074218751 -20.979392255760608,-69.02160644531251 -20.053351219365318)',4326)) as polygon,datasources.areas_desarrollo_indigena_2017_03_07_15 as datasrc) as sharedPolygon WHERE ST_IsEmpty(sharedPolygon)='false'"

Notice I've just prepended this: docker run "geodata/gdal:2.1.3" (I've also tried specifying the image without quotes, so that's not and issue).

This command complains it cannot connect to DB and also (see last error message line) interprets my Postgresql connection string as another datasource.

This is the output when I run the command through Docker:

ERROR 1: PQconnectdb failed.
could not connect to server: Connection refused
    Is the server running on host "localhost" (127.0.0.1) and accepting
    TCP/IP connections on port 5432?
could not connect to server: Cannot assign requested address
    Is the server running on host "localhost" (::1) and accepting
    TCP/IP connections on port 5432?

FAILURE:
Unable to open datasource `PG:host=localhost port=5432 dbname=mapserver user=postgres password=' with the following drivers.

Also it doesn't honor --debug ON parameter, it just spits previous string followed by the formats it supports.

UPDATE: Just found that I need to connect local database with container through Docker parameter --net, so the new command is:

docker run --rm --net=host -v /tmp/1:/data "geodata/gdal:2.1.3" ogr2ogr --debug ON -f 'KML' /tmp/1/areas_desarrollo_indigena_2017_03_07.kml PG:"host=localhost port=5432 dbname=mapserver user=postgres password=" -sql "SELECT * from ( SELECT ST_Intersection(datasrc.geom, polygon) as sharedPolygon, datasrc.* from ST_MakePolygon(ST_GeomFromText('LINESTRING(-69.02160644531251 -20.053351219365318,-69.78515625000001 -20.45789615504684,-69.09851074218751 -20.979392255760608,-69.02160644531251 -20.053351219365318)',4326)) as polygon,datasources.areas_desarrollo_indigena_2017_03_07_15 as datasrc) as sharedPolygon WHERE ST_IsEmpty(sharedPolygon)='false'"

However now it cannot write the output file despite the -v /tmp/1:/data, now it says:

PG: DBName="mapserver"
PG: PostgreSQL version string : 'PostgreSQL 10.1 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 7.2.0, 64-bit'
PG: PostGIS version string : '2.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1'
GDAL: GDALOpen(PG:host=localhost port=5432 dbname=mapserver user=postgres password=X this=0x23d1a20) succeeds as PostgreSQL.
GDAL: GDALDriver::Create(KML,/tmp/1/areas_desarrollo_indigena_2017_03_07.kml,0,0,0,Unknown,(nil))
KML: Attempt to create: /tmp/1/areas_desarrollo_indigena_2017_03_07.kml
ERROR 4: Failed to create KML file /tmp/1/areas_desarrollo_indigena_2017_03_07.kml.
ERROR 1: KML driver failed to create /tmp/1/areas_desarrollo_indigena_2017_03_07.kml
GDAL: GDALClose(PG:host=localhost port=5432 dbname=mapserver user=postgres password=X this=0x23d1a20)


texts:
I have a kml file with boundary polygon. It has a lot of gaps and overlaps in them. I converted them into shapefile and tried to use smooth and generalization feature tool in a couple of GIS software's like ArcGIS Desktop and Alteryx and didn't get any difference.

How could I smooth it up?


texts: Why are there holes in the resulting map of median reflectance values from Sentinel2 imagery on GEE?

texts:
I don not understand why are there holes in the resulting map of median reflectance values from S2 imagery on GEE and how do I fix it?

Also, why are all the clouds not getting masked?

//Function for masking clouds
// http://xzsunbest.tk/2018/07/05/MaskingOutCloudInSentinel2WithGEE/
// Bits 10 and 11 are clouds and cirrus, respectively.

var cloudBitMask = ee.Number(2).pow(10).int();
var cirrusBitMask = ee.Number(2).pow(11).int();
// Simple cloudMask function for Sentinel-2 based on QA band
function cloudmask(image) {
var qa = image.select('QA60');
// Both flags should be set to zero, indicating clear conditions.
 var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
       qa.bitwiseAnd(cirrusBitMask).eq(0));
       return image.updateMask(mask);
}

var filtered =  sentinel2.filterBounds(asset).filterDate('2015-01-01','2015-12-31'); //asset is Papua and West Papua i.e Indonesian side of PNG
print('Filtered (by date and ROI) Image Count: ', filtered.size()); //247

// apply the function on the image collection
var clipped = filtered.map(cloudmask);
print('Filtered and Masked Image Count: ', clipped.size());//247

// Make a composite image
var medianComp = clipped.median();
Map.addLayer(medianComp, {bands:['B4','B3','B2'], min: 0, max: 3500}, 'true color median composite', 0);


texts: Randomly filtering points in unevenly distributed point dataset?

texts:
I have a vegetation plots dataset which has unduly high concentrations of plots in certain parts of the dataset. I want to even out the concentrations of points so that there is a more even distribution of data geographically.

Can you point me to any previous questions that might have asked this question already or show me a process how I can use a grid and a point density map to filter the points to the desired maximum vegetation plots density (i.e. so many vegetation plots per sq km)?


texts: Concentric Ellipses around polygons as a representation of density

texts:
I have someone who has given me a thumbnail of a map with some concentric ellipsoids placed in estimated high density locations for our business. One location has three, and two other location have two and one respectively.

The assumption is that within the first ellipsoid is about 75% of our business, and within the second and third combined make up about 20% percent of the remaining.

They've asked me to create this, with numbers to back it up. The data I have right now is Zip Code Tabulation Areas for our geographical region, with registration information in a wide range of categories.

I've recently been utilizing GIS and delving into the world of spatial analysis, so forgive me if I'm naive. My thought was to find a technique that would allow me to fit an ellipsoid polygon around the polygons or their centroids for the density of ZCTA's that represent the 75% respectively, and create another one encompassing it to capture the remaining 20%. Those would be the areas we would focus on most heavily.

The question is... how would I even approach this in an application like QGIS? Furthermore, how would I break out the problem so that the two separate regions make up some of the calculations of that 75% and 20%.

Is this just something I have to approximate? My first thought would be to utilize some sort of binning technique like hexagon binning, but I think they're stuck in their desire to visualize it in terms of ellipses.

EDIT:

So I have a non-ideal solution I did in the absence of a proven, more automated technique.

I created a Points layer to create my ellipses, and used layer styles to created ellipses that fit roughly to the areas they indicated on their thumbnail, and with calculations fit the percentages they had estimated. It was all just a series of layer styles and manual selections from there to provide data with the summaries of coverage.

Not a great way in my opinion, but it gets the job done I suppose. I also learned that QGIS has terrible support for creating ellipses, even with plugins.


texts: How should I write the TimeBand for making a chart of accumulated precipitation between dates with GSMaP

texts:
I am trying to generate in GEE a chart of the accumulated rainfall between dates (two months approximately) for a series of polygons that I have in a shape file.
I am using the information from the "hourlyPrecipRate" band of the GSMAP satellite, but I am having problems with the time band.
How can I do this?
//Rain
var lluvia = ee.ImageCollection('JAXA/GPM_L3/GSMaP/v6/operational')
.map(function(im){return im.clip(Uy);}).filter(ee.Filter.date('2020-08-01', '2020-09-05'));
print(lluvia.size());
var precipitation = lluvia.select('hourlyPrecipRate');
var viz1 = {min: 0.0, max: 500.0, palette:['1621a2', 'ffffff', '03ffff', '13ff03', 'efff00', 'ffb103', 'ff2300'],};
var cumPrec = precipitation.sum();
Map.addLayer(cumPrec, viz1, 'cum Prec');

//Making the timeband for charts, here i need help
var withTimeBand = precipitation.map(function(image){
var datehours = ee.Image(ee.Number(image.get('system:time_start'))).toInt64();
return image.addBands(datehours.rename('datehours'));});
var minDates = dateHours.select('dateHours').min();
var maxDates = dateHours.select('dateHours').max();
var diff = maxDates.subtract(minDates);

//Should i get the data in days?
var diffDays = diff.divide(24);

//Making a chart with accumulated precipitation for all polygons
print('All_polygons',ui.Chart.image.series(precipitation, cz, ee.Reducer.first(), 500));

//Making a chart with accumulated precipitation for a specific polygon
print('Sp_polygons',ui.Chart.image.series(precipitation,cz.filterMetadata("Nombre","contains","Carlitos") , ee.Reducer.first(), 500));


texts: Can CartoDB read data that are hosted locally?

texts: Can CartoDB read from tables/data hosted on a local machine and not necessarily uploaded to  CartoDB's website?


texts:
I am having problem creating a query. I have a data map, which is organised in areas with codes (which in turn correspond to funding oportunities). I am trying to select for each code a layer. Up until then, its pretty straight forward.

The problem arises when I want to combine two codes. Through the Attribute Table, I select features using a expression,from the FIELDS AND VALUES, for example

"BINDUNG" LIKE '%11Z%'

However when I try to use

"BINDUNG" Like '%11Z%' AND LIKE '%812%'

or

"BINDUNG" Like '%11Z%' AND NOT LIKE '%812%'

it simply doesn´t work.

Is there another expression I can use?


texts:
I'm trying to georef old maps från 1960 that have clear boundary lines and my employer want an almost exact match with today's boundaries.

But it's very hard with 1th, 2th and 3th transformation.
With adjust och spline it's possible. But my boss says it will change the image and it can't be used.

Do spline or adjust change the original image?
How does they work?


texts: Dashes instead of "to" in mapinfo thematic legend

texts: The thematic legends for ranges are automatically set to the word "to", for example "10 to 20" "20 to 30" etc. I need to write it with a dash instead: "10-20", "20-30". Can I change the default settings somewhere?


texts: How to draw a raster over the World using Postgis

texts:
Im currently working on a system to track user interaction on mapwidgets. Now i need to visualize the bounds that users are looking at on the map. Since were using a whole lot of data i need to find a way of aggregating the data on the Server.
Im really just beginning to actually user the Postgis extension so im completely clueless how to do the following:
I thought about tiling the map into grids, starting at a few Hundred kilometers, then gradually tiling down those grids until i have grids of a few meters (like 110/70).
Now i thought id create a table:
ID (PK) | NE (Point) | SW (Point) | Diagonal (int or float) | Diagonal-1 (int or float)

So id save the northeastern Point the Southwestern Point, die Distance of Both of these Points (Diagonal) and the Distance of two ne/sw Points of the next smaller Raster. And do this for every Rectangle Starting at x:y then x/2:y/2 and so on until i have grids spanning the world of a few hundred meters. ( I guess creating a table will not serve me well since there will be a huuge amount of rows in the table. )
Now lets say the Server recieves an event of someone viewing the map.
bounds: ne,sw, center: point, diagonal: 3 kilometers

Id like to query for the ID of the Raster representing the bounds of the event.
so id check the grid table for the rows:
grid.diagonal >= events.diagonal > grid.diagonal-1

So an Event with a diagonal smaller then raster of size X but bigger then the next smaller raster will potentially get saved into a Raster of size X
then id like to check: event.center within which raster, and save the id of the raster in my event table. When querying the event table i could just do a count to get the number of events within that raster.
My guess is that I am looking at the problem in the wrong angle, has someone some ideas for me or done something similar like this?
edit, Well scratch the above lets make it a bit more simple
In essence i need the following:

I need a Grid spanning the World, and if custom viewbounds are supplied i need to check which Grid these viewbounds touch and save the bounds of that grid.


texts: How to set a map's extent to the envelope in miles

texts:
I want to set extent of the map into envelope which sizes are defined by user in miles(kilometers, whatever human readable unit). How to make it? I wrote the following code:

static void ZoomToPoint(IActiveView focusMap, ExportJson2Feature.Geometry geom)
{
    IPoint point = new PointClass();
    point.PutCoords(geom.X, geom.Y);

    IEnvelope newExtentEnvelope = point.Envelope;
    newExtentEnvelope.Expand(1, 1, false);
    ((IGeometry)newExtentEnvelope).Project( ((IMap)focusMap).SpatialReference );
    focusMap.Extent = newExtentEnvelope;
}

How to make it know that sizeX/sizeY are in miles?


texts: Create a new NetCDF file combining information from two other NetCDF files in R

texts:
I wish to combine the "t" variable from the second .nc file, and use the lat, lon variables from the first .nc file as the lat and lon of the t variable, and store them to a new .nc file, that will have as variable the "t" and as lat and lon the "lat" and "lon" variables from the second .nc file. Until now, this is what I have managed to do:

library(ncdf4)
coords.nc = nc_open("coordinates.nc")
lat = ncvar_get(coords.nc, "lat")
lon = ncvar_get(coords.nc, "lon")
nc_close(coords.nc)

temp.nc = nc_open("temp.nc")
t = ncvar_get(temp.nc, "t")
nc_close(temp.nc)

# Define the dimensions and variables
dimX = ncdim_def("Long", "degrees", lon)
dimY = ncdim_def("Lat", "degrees", lat)
dimT = ncdim_def("t", "mm", pr)
varT = ncvar_def("t", "mm", dimT, t, prec = "float")

# Define missing value
mv = -9999

# Define the data
var1d = ncvar_def( "var1d", "units", dimX, mv, prec="float")
var2d = ncvar_def( "var2d", "units", list(dimX,dimY), mv, prec="float")
var3d = ncvar_def( "var3d", "units", list(dimX,dimY,dimT), mv,   prec="float")

# Create the NetCDF file
my_new_nc = nc_create("writevals.nc",  list(var1d, var2d, var3d))

# Write data to the NetCDF file
ncvar_put(my_new_nc, var3d, varT, start=c(1, 1, 1), count=c(1, 1, length(varT)))

# Close the new file to finish writing
nc_close(my_new_nc)

The code works fine until I reach the line for "nc_create":

nc = nc_create("writevals.nc",  list(var1d, var2d, var3d))

the error message I get is the following:

Error in R_nc4_enddef: NetCDF: One or more variable sizes violate  format constraints
Error in R_nc4_sync: NetCDF: Operation not allowed in define mode
Error in Rsx_nc4_put_vara_double: NetCDF: Operation not allowed in define mode
Error in ncdim_create(nc, d, verbose) :
Error in ncdim_create, while writing dimvar values!

Could you please help me solve it?


texts: Copying content from one Portal for ArcGIS to another using Python?

texts: I am trying to copy content from one portal to another with python only, without using AGO assistant https://ago-assistant.esri.com/, as this doesn't work behind our firewall in our test network. I have explored these links for obtaining the correct python script to be able to perform the task: https://github.com/Esri/ago-tools/tree/master/samples but I don't know which script works for the copying of content from one to another? If anyone has used these, simply pass on the name of the script and that will suffice this request?


texts: Keep synchronization when simplify features from different data source

texts:
We have two shapefile data sources:

1 water-riverbank.shp

With geometry type of Polygon which stands for the area of water

2water-riverline.shp.

With geometry type of LineString which stands for the direction of the water. We use these features to label the water.

Now we use mapshaper to simplify these data, however we may get this kind of result:

The simplified riverbank and riverline are not synchronized.

Any way to avoid this?


texts: Getting directions using routes in ArcGIS API for JavaScript?

texts:
I am creating several routes, these routes represent the public transportation network of a city. Now I want to calculate the best combination to go from a point A to a point B using those routes.

Right now I am working with ArcGIS API for JavaScript 3.22, but I am open to use there others SDK from Arcgis.

I'm checking the docs, and the samples, but I have no clue about how to implement this. I am trying to use the RouteTask, but it looks like I'll need to use another component.


texts: Is there a basic hands-on example of the use of UK GIS data (not in the ESRI native format) being imported for use in ArcMap?

texts:
Is there a good, free, generally available hands-on UK-based example of UK GIS data (not in the ESRI native format) being imported for use in ArcMap, using ESRI tools only?

I want to find and organize GIS data to create a map that ArcGIS Desktop can consume.  The example would need to furnish some data and maps, rather than assuming that I already have something suitable.

I am not looking for a solution to a programming problem; I am hoping to demonstrate that ArcMap can work with, for example, OS data.


texts: Deleting data by location using another shapefile

texts:
I have two shapefiles: one is main file and another is a part divided (manually by selecting) from the main file.
Now i want to delete the divided part from the main file with the help of divided part. How this possible in QGIS?


texts: How display tiles by using earth plugin in browser?

texts:
I was using MapTiler to generate tiles for my maps, and everything is great. But I have problems with displaying those tiles in 3D version by using earth plugin. The data which gave me MapTiler was with files but for API 2, but I want to use now API 3 version.

All works fine with using maps, but I cannot display tiles by using earth plugin.

Please give me some examples, if anyone had to do something with this subject.


texts:
I am trying to generate subbasin with provided multiple outlets points by using r.stream.basin. I have install r.stream.basin addon in grass but it does not appear in qgis. Is there any ways to use this addon in qgis??

answer posted here link does not help.


texts: Using ST_split to split a line feature using nearby points

texts:
I have a line shapefile and a point shapefile, and

the line has just one segment

and its attribute table,

the point shapefile has points near the line above, not on it

and the close look,

I want to use ST_split to split this line at locations which are nearest to the points I want to get the original line with more segments(because being split). For example, find the location on the line, which is the nearest to the point, and split the line at this location(star).

I loaded this line shapefile into PostGIS and here's my SQL code in PostgreSQL,

select ST_AsText(ST_Split(ST_GeomFromText('LINESTRING(0 0, 247907.60 2746404.62, 200    200, 400 400)'),
       ST_Line_Interpolate_Point(ST_GeomFromText('LINESTRING(0 0, 247907.60 2746404.62, 200 200, 400 400)'), 0.2)))

And I checked the table of this line, it wasn't being splited.

Where did I make mistakes？ I am new to use ST_split, please give me any suggestion, thanks.

I am using QGIS 2.2, PostgreSQL 9.2 and PostGIS 2.0.4 under Windows 8.1 x64.

Update#1

@simplexio, I used your code,

select ST_Line_SubString(ST_GeomFromText('LINESTRING(0 0, 247907.60 2746404.62, 200 200, 400 400)'), 0 , 0.200 ) as geom into split_line ;

and got a new table named split_line,

But as I used the following code you posted,

SELECT ST_Line_SubString(line, 0, ST_Line_Locate_Point(l.line, p.point)) as from_start_geom,
       ST_Line_SubString(l.line,ST_Line_Locate_Point(l.line, p.point), 1) to_end_geom, l.line_id, p.point_id FROM lines l , points p
WHERE l.line_id = 1 , p.point_id = 1

I got a syntax error,

How to fix this and move on？ And where should I assign the line I want to split as input of st_split？


texts: Showing layers features data on click in OpenLayers5?

texts:
I am trying to display layer data on map click but I can't figure out how to do it. There is not an example of this important feature.

I believe it's the getfeatureinfo event. Here is the messy API documentation. Here is what i am trying i am getting js.js:148 Uncaught ReferenceError: view is not defined

Update:

I am following this example only and i managed to make it work.
Though, it's not what I actually want to achieve. I want to access columns as variables and not display the whole table in an iframe as shown in this poor example

This is closer of what i am trying but it loads an xml file. Another weird example.

    var map = new ol.Map({
      target: 'map',
      layers: [baseLayer, wmsLayer],
      view: new ol.View({
        center: ol.proj.fromLonLat([33.33386, 35.14710]),
        zoom: 14
      })
    });
       map.on('click', function (evt) {
  document.getElementById('info').innerHTML = '';
  var viewResolution = map.getView().getResolution();
  var url = wmsSource.getGetFeatureInfoUrl(
    evt.coordinate, viewResolution, 'EPSG:3857', {
      'INFO_FORMAT': 'text/html'
    });

  if (url) {
    console.log(url);
    document.getElementById('info').innerHTML = '<iframe seamless src="' + url + '"></iframe>';
    // fetch(url).then(function (response) {
    //   return response.text();
    // }).then(function (response) {
    //   // this is the standard way to read the features
    //   var allFeatures = new WMSGetFeatureInfo().readFeatures(response);
    //   console.log(allFeatures.length.toString());
    // });
  }
});


texts: What WGS84 epoch date is used by Google Maps

texts:
I have some Nad83 data that I want to overlay on Google Maps.  This article by ngs.noaa.gov explains translating from nad83 to wgs84 better than anything I have read. If I were to use a tool like the noaa transformation tool to convert a coordinate pair from NAD83 to WGS84 for this purpose, would I want to use an output epoch that matches the approximate date of the imagery?  Or is all google imagery georeferenced based on some common epoch independent of the date the imagery was created?

I realize that this datum shift will be small and probably hard to realize on the +/- 1 meter Google photography.  But still, I want to make sure I get the best alignment possible.


texts: Recommended ArcGIS Server and RDBMS layout on these two machines?

texts:
We're looking at setting up ArcGIS Server with Postgres for the RDBMS. We have two machines, one aging and a new hotness. What is the best implementation strategy with these resources? Should we put both Server and Postgres on the new machine, put both on both machines, Server on new and Postgres on old, Postgres on new Server on old? ...or something else entirely?

Old and not yet busted:

Year: 2008
CPU: 2x Dual-core Opteron 2.8GHz
RAM: 64GB
Network: 2x 1Gbps
Storage: 1.3TB in RAID 5 (6x 280GB 15K rpm serial SCSI)
OS: Windows Server 2003 x64 (upgrade to Server 2012 planned)

New hotness:

Year: 2014
CPU: 1x 8-core @2Ghz (Xeon E5-2640 v2)
RAM: 16GB (2x 8GB 1600Mhz DDR3L)
Network: 2x 1Gbps
Storage: 3.6TB in RAID 5+0 (6x 900GB 10k rpm serial scsi)*
OS: Windows Server 2012 x64

Storage: this is projected and open to change, we haven't purchased the drives yet.

File resources which will be put into the RDBMS. Expected data growth rate: 15%/yr (but really, we have no idea):

File-gdb:   15gb
Rasters:  1000gb (includes overviews)

Internal users, same subnet, mostly gigabit:

6 to 10 active daily use ArcGIS Desktop users drawing from the data,
mostly read operations, pulling from the central repository for base
maps and shared data, adding their own stuff from elsewhere.
2 to 3 data editors, but this happens rarely, once a week perhaps.
25-50 readers, currently using ArcReader but looking to migrate to
web maps (a major reason we're implementing ArcGIS Server to begin
with).

Projected purpose of ArcGIS Server + RDBMS:

primary daily use warehouse for the users mentioned above
Master server for geodatabase replication to district offices
development and learning area for web maps and services which will eventually be migrated to outside firewall for public consumption (which is ArcGIS Server + Oracle)
(ArcGIS Portal?)

We have whatever licenses for ArcGIS Server needed (from Basic to Enterprise Advanced), but not for Oracle or SQL Server, hence the Postgres for RDBMS specification.

There is no Sys Admin IT or DBA staff to dedicate to setting up and maintaining these machines and services, that will be handled by a GIS person who still needs to allocate most of his/her time to regular ol' GIS. (guess who?) This is a preface to saying the solution that makes best use of the hard and software is not necessarily the optimal solution, please keep administrative overhead in mind.


texts:
Is it possible to generate map tiles from QGIS?

I guess QGIS Server must have a way to do it, but is there a way to just generate the tiles and save them to disk?

With the new symbology renderer and scale-dependent rendering, QGIS is a powerful tool for generating dynamic maps.

Since I'm using the new-symbology renderer, qgis-mapnik is not an option, and since I'm using scale-dependent rendering, just generating a large georaster and using gdal2tiles to tile it is also not an option.


texts:
How can I find features that do not touch any other feature? I have about 4000 pipeline elements and I want to know which of them are separated from the network. I tried with:

SELECT h1.gid FROM pipelines h1, pipelines h2 WHERE ST_Disjoint(h1.geom, h2.geom)

but I get too many results.


texts: Using Leaflet, PostGIS and pgRouting for web mapping

texts:
I am new to Open Source GIS. I have already created a web mapping using custom tiles and Leaflet. I have used my own database. TileMill was used to create the tiles and Leaflet was used to visualize the tiles. Now I want to create a query of routing and Drive direction on the web mapping. I saw some tutorials of underdark in the following blog where routing visualization was done through QGIS software.
http://underdark.wordpress.com/tag/pgrouting/
Now I want to know that whether it is possible to use my custom data (shape file) as PostGIS to visualize routing on web map or not?


texts: Select layer by attribute for entries in a dictionary (python)

texts:
I'm kinda new to python (and programming) and having trouble figuring out how to properly specify what I want to do.

Basically, my code compares names of well owners and orps owners and saves all orps that match a given well in a dictionary as

matches[fid of well] = fid of matching orps, fid of matching orps, etc.

Such that some wells are not included in the dictionary, some wells have one entry, and some wells have several entries.

Once created, I then need to:

select the well that the key of the dictionary is referring to
select the orps pts that the key is referring to
create a near table of matching orps for that given well. (this part I know how to do)

I guess I'm just unsure of the syntax required to do this. I'm sure the answer to this lies somewhere within the SQL query in the:

arcpy.SelectLayerByAttribute_management(in_layer, {selection_type}, {where_clause})

function, but I'm having trouble constructing the query. Can anyone give me a hand?

Full code:

import arcpy

arcpy.env.overwriteOutput = 1

# define a workspace
arcpy.env.workspace = r"C:\Users\tmc18\Desktop\comp_orps\NYS_Wells.gdb"

# Define input files
orps = r"C:\GIS\Madison_County_orps.gdb\orps_Cazenovia"
wells = r"C:\GIS\NYS_Wells.gdb\Madison_wells_cazenovia"

# Make a feature layer
arcpy.MakeFeatureLayer_management(orps, "orps1")
arcpy.MakeFeatureLayer_management(wells, "wells1")
wellsFL = "wells1"
orpsFL = "orps1"

# Create dictionary of last names of all well owners
well_owners = {}
rows = arcpy.da.SearchCursor(wellsFL, ["OBJECTID", "OWNER_L"])
for row in rows:
    well_owners[row[0]] = str(row[1])
del rows

# create dictionary of addresses of ORPS pts
loc_info = {}
rows = arcpy.da.SearchCursor(orpsFL, ["OBJECTID", "LOC_NUM", "LOC_NAME", "PRINT_KEY"])
for row in rows:
    loc_info[row[0]] = str(row[1]+" "+row[2]+","+" "+row[3])
del rows

matches = {}

y = 0

# Create search cursor which will iterate through wells
with arcpy.da.SearchCursor(wellsFL, ["OBJECTID"]) as cursor:
        for row in cursor:

            match_list = []
            # set well owner name for this record
            record = row[0]
            owner = well_owners[record]

            # create search cursor which will iterate through orps owners
            with arcpy.da.SearchCursor(orpsFL, ["OWNER_L", "OBJECTID"]) as orpscur:
                for row in orpscur:
                    # set orps owner name
                    orpsowner = row[0]
                    ID = row[1]

                    # compare owner's names
                    if owner != orpsowner:
                        pass

                    else:
                        y = y + 1
                        match_list.append(ID)

            matches[record] = match_list

print y, "total matches found"


texts: How to publish a layer with standalone geowebchache in geoserver?

texts: I have a multiple geoservers on two different systems with same port(one is backup for another) and geowebcache in on of the system with different port . Now i want to publish my layers in  geoserver with  standalone geowebcache. is there any material to refer? Please help me?


texts: Kernel Density Estimate on point data in degrees (scikit-learn)

texts:
I am estimating point density, where I am having point coordinates in degrees, I need a raster(over the extent of these points) which provides a density estimation on each cell.

I have tried on R but I have not got any success as the extent of my point data is very large. I am trying scikit-learn package in python, but stucked with the implementation.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity
import matplotlib.pyplot as plt
import gdal
import os
import ogr
import osr
from sklearn.cluster import DBSCAN

def density_Viz(str_InputSHP):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    ds = driver.Open(str_InputSHP)
    layer = ds.GetLayer(0)
    lst_train = []
    lst_coord = []
    for feat in layer:
        geom = feat.GetGeometryRef()
        X = geom.GetX()
        Y = geom.GetY()
        xy = (X, Y)
        lst_coord.append(xy)
        lst_train.append([X, Y])

    # val=getKernelDensityEstimation(lst_coord)
    lst_train = np.array(lst_train)
    lat_lon_val = lst_train.copy()
    lst_train *= np.pi / 180.
    # create grid
    (dblXmin, dblXmax, dblYmin, dblYmax) = layer.GetExtent()
    pixelWidth = 0.00082
    Nx = int((dblXmax - dblXmin)/pixelWidth)
    Ny = int((dblYmax - dblYmin)/pixelWidth)
    xmin = dblXmin + 0.0008
    xmax = xmin + (Nx * pixelWidth)
    ymin = dblYmin + pixelWidth
    ymax = ymin + (Ny * pixelWidth)
    # x coordinates of the grid cells
    xgrid = np.arange(xmin, xmax, pixelWidth)
    # y coordinates of the grid cells
    ygrid = np.arange(ymin, ymax, pixelWidth)

    X, Y = np.meshgrid(xgrid, ygrid)
    xy = np.vstack([X.ravel(), Y.ravel()]).T
    xy *= np.pi / 180.

    kde = KernelDensity(bandwidth=(0.0032*np.pi/180), metric='haversine',
                        kernel='epanechnikov', algorithm='ball_tree', breadth_first=False).fit(lst_train)

    val = kde.score_samples(lst_train)
    val1 = kde.score_samples(xy)

    print(np.amax(val))  # getting 16, expected value is nearly 70
    print(np.amin(val))  # getting 13
    print(np.amax(val1))  # getting 16
    print(np.amin(val1))  # getting -Inf
    print("Done")


texts: GIS Action that loads a raster, creates pyramids and sets color

texts:
Using the code in "Developing QGIS Action which loads a raster".

from PyQt4.QtCore import QFileInfo,QSettings
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem

s = QSettings()
oldValidation = s.value( "/Projections/defaultBehaviour" )
s.setValue( "/Projections/defaultBehaviour", "useGlobal" )
vl = QgsMapLayerRegistry.instance().mapLayersByName('catalogo_CTR')[0]
qgis.utils.iface.setActiveLayer(vl)
root = QgsProject.instance().layerTreeRoot()
group_name = "Raster layers"
group = root.findGroup(group_name)
if group == None:
    group = root.addGroup("Raster layers")
else:
    pass
fileName = '\\\\rvphnas02pw\Prodotti_Cartografici\\Ortofoto\\2012_AGEA\\TIFF_GBO\\[% A_CODICE %].tif'
fileInfo = QFileInfo(fileName)
baseName = '[% A_CODICE %]'
rlayer = QgsRasterLayer(fileName, baseName)
crs = QgsCoordinateReferenceSystem()
crs.createFromSrid(3003)
rlayer.setCrs(crs)
QgsMapLayerRegistry.instance().addMapLayer(rlayer, False)
group.insertChildNode(-1, QgsLayerTreeLayer(rlayer))

s.setValue( "/Projections/defaultBehaviour", oldValidation )

Is it possible to create automatically the image-pyramids and set color white as trasparent?


texts: SLD one layer with several feature types

texts:
I wrote an SLD for a layer that contains two feature types: us-net-sw:SewerPipeand us-net-common:Appurtenance. My WMS - implemented with deegree - only renders us-net-sw:SewerPipe. I don't know if my sld is correct, I'm trying to define one layer with two FeatureTypeStyle. If I remove the us-net-sw:SewerPipe FeatureTypeStyle deegree renders the us-net-common:Appurtenance.

<?xml version="1.0" encoding="utf-8"?>
<sld:StyledLayerDescriptor xmlns:sld="http://www.opengis.net/sld" xmlns:app="http://www.deegree.org/app" xmlns:us-net-sw="http://inspire.ec.europa.eu/schemas/us-net-sw/4.0" xmlns:deegreeogc="http://www.deegree.org/ogc" xmlns:ogc="http://www.opengis.net/ogc" xmlns:us-net-common="http://inspire.ec.europa.eu/schemas/us-net-common/4.0" xmlns:sed="http://www.deegree.org/se" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:se="http://www.opengis.net/se" xmlns:net="http://inspire.ec.europa.eu/schemas/net/4.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd">
    <sld:NamedLayer>
        <se:Name>US.SewerNetwork</se:Name>
        <sld:UserStyle>
            <sld:IsDefault>1</sld:IsDefault>
            <se:FeatureTypeStyle>
                <se:FeatureTypeName>us-net-sw:SewerPipe</se:FeatureTypeName>
                <se:Rule>
                    <se:LineSymbolizer>
                        <se:Geometry>
                            <ogc:PropertyName>net:link/us-net-common:UtilityLink/net:centrelineGeometry </ogc:PropertyName>
                        </se:Geometry>
                        <se:Stroke>
                            <se:SvgParameter name="stroke">#000000</se:SvgParameter>
                            <se:SvgParameter name="stroke-width">1</se:SvgParameter>
                        </se:Stroke>
                    </se:LineSymbolizer>
                </se:Rule>
            </se:FeatureTypeStyle>
            <se:FeatureTypeStyle>
                <se:FeatureTypeName>us-net-common:Appurtenance</se:FeatureTypeName>
                <se:Rule>
                    <se:PointSymbolizer>
                        <se:Geometry>
                            <ogc:PropertyName>net:geometry</ogc:PropertyName>
                        </se:Geometry>
                        <se:Graphic>
                            <se:Mark>
                                <se:WellKnownName>square</se:WellKnownName>
                                <se:Fill>
                                    <se:SvgParameter name="fill">#808080</se:SvgParameter>
                                </se:Fill>
                                <se:Stroke>
                                    <se:SvgParameter name="stroke">#000000</se:SvgParameter>
                                    <se:SvgParameter name="stroke-width">1</se:SvgParameter>
                                </se:Stroke>
                            </se:Mark>
                            <se:Size>5</se:Size>
                        </se:Graphic>
                    </se:PointSymbolizer>
                </se:Rule>
            </se:FeatureTypeStyle>
        </sld:UserStyle>
    </sld:NamedLayer>
</sld:StyledLayerDescriptor>


texts: Identifying buildings by address in OpenStreetMap with QGIS?

texts:
I'm trying to find a way to represent San Francisco's list of softstory buildings on OSM as building polygons. The city of San Francisco produces a list of addresses as a Fusion table. The points produced in Fusion are not what I need though.

I am looking for a way to build a thematic map of the buildings that need to be retrofitted. I'm playing with QGis, importing the building features using TurboPass but I'm not sure which features I need to import. At first I tried using MMQGIS geolocation thinking that I could just select building's polygons that included the points but many of the points end up in the middle of a road or far from the actual building.

What OpenStreetMap tags and keys should I query and import in QGis to get the full addresses of buildings?


texts: I tried to use QGIS's Distance to Nearest Hub to calculate the distance from observation points to the nearest stream (I had converted my stream line layer to a point layer with points close together, as suggested in another post here). I chose meters as my units of the output file. The resulting table had a field called HubDist, but the values were all 0.00XXX, such as 0.00415, 0.00513, etc, which didn't make sense if the units were truly meters as the observation points should have been typically in the 100-1,000 meter range. However, when I used the map ruler to measure some of the points to the nearest streams, it appears as if this might be a decimal point issue. That is, the 0.00415 point seemed to be about 415 meters away from the stream. Maybe this is a strange coincidence. In any case, can someone please tell me what might be the cause of these distances being incorrect and how I might fix them? Right now, these two source layers are unprojected WGS 84; do I need to project them first? Or am I correct that I need to only move the decimal point to get units in meters, which seems strange.


texts:
When I run select pgr_createTopology('public.roads', 0.0001, 'geom', 'id');,
all I get is "FAIL". The roads table has these fields:
id (integer)
source (integer)
target (integer)
geom (geometry LineStringM,3006)
What am I missing? I have been unable to find a list of requirements for the network table, which I create using FME. I use pgRouting 2.2.2.


texts:
I would like to use marker symbols for markers on my Mapbox (v3.0.1) map, I have found several sources that say I can use icons from this page. Yes, some of them work for me (danger, religious-christian), but some of them do not, for example gaming or toilet.

Does anybody know the real available list of marker-symbols?


texts:
I'm having a problem with visualizing a WFS layer in Openlayers. I can see I get the response, but nothing shows up on my map. I tried more potential solutions from stackexchange, but couldn't get it working. The WFS layers are shown correctly in QGIS without any special configuration.
The address of the WFS in question is: http://eagri.cz/public/app/wms/plpis_wfs.fcgi - open data for the Czech Republic.
My code (basic version, I tried more solutions which didn't work):
var vectorSource = new ol.source.Vector({
    format: new ol.format.WFS(),
    url: function (extent) {
        return (
            'http://eagri.cz/public/app/wms/plpis_wfs.fcgi?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=LPIS_DPB_UCINNE&SRSNAME=EPSG:102067&' +
            'bbox=' +
            extent.join(',') +
            ',EPSG:3857'
        );
    },
    strategy: ol.loadingstrategy.bbox
});

var vector = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: '#FF0000',
            width: 2
        }),
        projection: new ol.proj.Projection('EPSG:102067')
    }),
});

var map = new ol.Map({
    target: 'map',
    pixelRatio: 1,
    layers: [
        vector
    ],
    view: new ol.View({
        center: ol.proj.fromLonLat(cw_global_initCoordinates),
        zoom: 15
    })
});

Generated request example:
http://eagri.cz/public/app/wms/plpis_wfs.fcgi?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=LPIS_DPB_UCINNE&SRSNAME=EPSG:102067&bbox=1885304.656759394,6323618.39452889,1892413.3003899155,6326422.678004101,EPSG:3857 - follow the request link to get the example response


texts:
I've problems symbolizing the result of a GeoJSON. I do the function and define the intervals based on the results of the GeoJSON. I want to symbolize the abundance field.
I put a part of the GeoJSON, example:
{"type":"Feature","geometry":{"type":"MultiPolygon","coordinates":[[[[-9,43.5],[-9,44],[-8,44],[-8,43.5],[-9,43.5]]]]},"properties":{"id":"3446","aphia_id_code":"140600","camp":"N13","abundance":"39.898","geojson":"{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-9,43.5],[-9,44],[-8,44],[-8,43.5],[-9,43.5]]]]}"},"id":18},{"type":"Feature","geometry":{"type":"MultiPolygon","coordinates":[[[[-5,43.5],[-5,44],[-4,44],[-4,43.5],[-5,43.5]]]]},"properties":{"id":"3838","aphia_id_code":"140600","camp":"N13","abundance":"2.205","geojson":"{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-5,43.5],[-5,44],[-4,44],[-4,43.5],[-5,43.5]]]]}"},"id":19}

And this is the styling function
var style = new ol.style.Style({
    stroke: new ol.style.Stroke({color: '#425364', width: 1}),
    fill: new ol.style.Fill()
});

var dist = function(feature) {
    var value = feature.get('abundance');
    var color = value <  100 ? '#F6CECE' :
                value <  1.500 ? '#F5A9A9' :
                value <  12.000 ? '#F78181' :
                value <  22.000 ? '#FA5858' :
                value < 43.000 ? '#FE2E2E' : '#ff0000';
    style.getFill().setColor(color);
    return style;
};

When I represent it, it does not respect the color code, it sets the same color for all intervals. I think if this can be caused by formatting numbers, but I'm not really sure.  How can I represent these intervals with those values?


texts: Converting KML to feature class using ModelBuilder?

texts:
I am experiencing issues tyring to convert a KML to a layer file to use within ModelBuilder. The KML contains a set of polygon data and a set of polyline data, the polygon data is the data I wish to work with. When I set up the "KML to Layer" within the Model I can't link it to the convert "Feature Class to Feature Class" tool and proceed with my model, it gets stuck at this point. What am I doing wrong?

I am using ArcMap 10.1


texts: Aerial images coordinate conversion problem from ArcMap to QGIS

texts:
I have multiple different historical aerial images from an ArcMap project and I want to "extract" somehow the present-day satellite image of their areas from Google Satellite.

So what I tried to do is, import Google Satellite with the XYZ Tiles Plugin, set the Coordinate Reference System (CRS) to DHDN / 3-degree Gauss-Kruger zone 4 (EPSG:31468), which is the CRS from the images, and saved the area with the Layer --> Save As functionality, using the coordinates of the images as north, east, south, west.

However, this does not work. I either get a blank (white) image or the image is of a completely different location and if I just import the image into QGIS it is also not at the right location.

Yet, when I import the image into ArcMap it is at the correct position.

The output of gdalinfo on one of the images is:

gdalinfo 1978-03.tif

Driver: GTiff/GeoTIFF
Files: 1978-03.tif
       1978-03.tif.ovr
       1978-03.tif.aux.xml
Size is 10766, 10619
Coordinate System is `'
Origin = (0.000000000000000,147.486111111111114)
Pixel Size = (0.013888888888889,-0.013888888888889)
GCP Projection =
PROJCS["DHDN_3_Degree_Gauss_Zone_4",
    GEOGCS["GCS_Deutsches_Hauptdreiecksnetz",
        DATUM["Deutsches_Hauptdreiecksnetz",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",4500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",12.0],
    PARAMETER["Scale_Factor",1.0],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0],
    AUTHORITY["EPSG","31468"]]
GCP[  0]: Id=, Info=
          (8924.44718310842,9456.67887323909) -> (4469283.68817037,5535422.30978963,0)
GCP[  1]: Id=, Info=
          (9056.91760563757,9337.02816901542) -> (4469314.11531456,5535445.32858567,0)
GCP[  2]: Id=, Info=
          (7505.95437375864,7683.60903820582) -> (4469025.98348829,5535792.51553005,0)
GCP[  3]: Id=, Info=
          (7620.95041533496,7634.67871375615) -> (4469048.63186693,5535802.46388328,0)
GCP[  4]: Id=, Info=
          (10292.3337961421,8017.17334593507) -> (4469579.09605285,5535697.31825632,0)
GCP[  5]: Id=, Info=
          (10244.4309673187,8520.11818380561) -> (4469562.48018629,5535595.82388666,0)
GCP[  6]: Id=, Info=
          (6920.10031108643,7016.7363824416) -> (4468912.66221999,5535935.89352514,0)
GCP[  7]: Id=, Info=
          (1976.54547536775,749.573714908212) -> (4467982.38535943,5537290.12192803,0)
GCP[  8]: Id=, Info=
          (6453.23938492237,865.984278769698) -> (4468902.69622784,5537184.68525883,0)
Metadata:
  TIFFTAG_DATETIME=2015:05:11 07:26:38
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS4 Windows
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (       0.000,     147.486)
Lower Left  (   0.0000000,   0.0000000)
Upper Right (     149.528,     147.486)
Lower Right ( 149.5277778,   0.0000000)
Center      (  74.7638889,  73.7430556)
Band 1 Block=10766x24 Type=Byte, ColorInterp=Gray
  Overviews: 5383x5310, 2692x2655, 1346x1328, 673x664, 337x332, 169x166

As I have zero experience with ArcMap, I thought about using QGIS for this task. Furthermore, I can only use ArcMap at my University.


texts: Merging small lines to form longer line in ArcGIS Desktop?

texts:
I'm having some trouble with the dissolve tool in arcgis. I have a stream network that has some streams made up of a number of little segment that I would  like to transform so  that all the stream portions that are connected to each other appear as one entity, a longer stream.

Here is a zoom with selected portions of streams in blue to illustrate my problem

I'd like to have all these portions of the stream (blue and black) as one entity but I also need to keep the informations such as length (long_m) for the whole stream once they are merged together. I assume that the statistics field with the statistic type SUM can make this work.

Here is what my attribute table looks like.

As I said I have tried the dissolve tool but haven't been able to find a way to get the streams like I wanted. Ideally each small part composing the stream would be added to each other considered as one entity and the length of each portions would be added to form one stream with one longer length.

here is how I set up my dissolve, I've also tried other to check unsplit lines but it didn't work.


texts: Remove first n characters in label using ArcGIS 10 label expression builder with VBscript

texts:
I have a planning cadastre dataset. It has the type of plan (DP, RP or SP) and then the plan number, thusly:

DP123456
RP1248590000
SP3855
etc...

I want to make a label in the expression builder in ArcGIS 10 which skips the first two characters and only has the plan number. In full, this is what I have done:

left([planlabel],2)&vbnewline&left([planlabel],len([planlabel]-2))&vbnewline& [lotnumber]

You can see here that I am making a stacked label which will put the type of plan on the first row (it is always the first two characters of [planlabel] and then the second row has the plan number, followed by [lotnumber] on the last line.

In Excel, string indexing of this sort is apparently achieved by CELL,len(CELL-2) and I was hoping the same would hold for VBscript in ArcGIS but apparently this is not the case.

Any clues?


texts: Extracting points as CSV from count polygon

texts:
By using "points count polygon" I found the number of points in every polygon. I have to extract those points to a CSV file but I can't.

I came to know that I have 957 points in my 1st polygon (as shown in figure attached) now I need those 957 points to be extracted as CSV file.

How can I do this?


texts:
I've had a look at the following thread about using MXD2QGS and so far I've figured out that this only works with version 10.0 and Python 2.6. I installed MXD2QGS on a machine with Arc10.0 and Python 2.6 and ran the tool.

It's throwing an error:

  No module named ext

I opened the python script and it imports xml.dom.ext, so I think this is where the error is occurring. I did some research and it seems to be that this module doesn't come packaged with 2.6 and I think I read that the last version it came with was 2.4. That's fine, but I'm wondering how the developer of this tool created it in 2.6? I couldn't find this module to download, so I'm also wondering if there's a workaround module?

If anyone has some suggestions, I'm all ears.

Thanks,
Mike


texts: gdal_array.SaveArray() leaves dataset open in Python

texts:
When I use gdal_array.SaveArray() to create a raster, the newly created dataset appears to stay open in Python, preventing other processes from working with it. For instance, consider the following (super minimal) code:

>>> a = np.arange(300).reshape((3, 10, 10))
>>> gdal_array.SaveArray(a, "test.tif")
<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x0458B968> >
>>> os.rename("test.tif", "test2.tif")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
WindowsError: [Error 32] The process cannot access the file because it is being used by another process

I am similarly prevented from moving, renaming, or deleting the file directly from Windows explorer, with the message The action can't be completed because the file is open in python.exe. More importantly, I can't open the file some visualization/processing programs that I use. Once I exit python, the file is "released" and I can manipulate it to my heart's content.

It doesn't seem like the file has a name associated with it, so I can't close it as I would a raster that I had specifically opened:

>>> dir()
['__builtins__', '__doc__', '__name__', '__package__', 'a', 'gdal', 'gdal_array', 'np', 'os']

What causes this behavior? Is there a way to call SaveArray() such that it doesn't keep the file opened after writing it? Or, a way to close the file from within Python?

In case it's important, my Python bindings are from gdal 1.11.1 for Python 2.7.8 on Windows 7.


texts: gxp_wmsgetfeatureinfo GetFeatureInfo on a non-advertised layer possible?

texts:
So, as the title states, I will like to know if it is possible to use gxp_wmscsource to send an GetFeatureInfo query to a non-advertised layer that was added to my opengeo SDK app using lazy loading.

I made my app and added a gxp_wmsgetfeatureinfo to it, but if a click over a non-advertised layer, the app doesn’t even sends a GET, I suppose it is because there is not a list of available layers...

I know it is possible to send a WMSGetFeatureInfo to a layer even if it is not advertised but I don’t know how to tell gxp_wmsgetfeatureinfo to do it.


texts: Projecting Irish TM65 to ITM95 in ArcGIS 10.6

texts:
Some years ago I created a polygon coverage of Wexford Harbour in Ireland by tracing a rectified aerial image of the harbour using ArcGIS. Four tic points on the imagery were known locations; their positions were converted to Irish national grid (TM65) using an early version of GridInQuest. By georeferencing the image, the locations in the traced shapefile corresponded to Irish Grid.

In order to use ArcGIS online mapping features I needed to use IRENET95 so using  its helpful guide I installed the package

  ntv2-grids-master IG to ITM.

Using ArcToolbox I defined the projection as TM65_IG and attempted to project to IRENET95_Irish_Transverse_Mercator. Everything appeared to go perfectly; the box below Geographic Transformation displayed:

  TM65_To_ETRS_1989NTv2 + IRENET95_To_ETRS_1989

When I pressed OK the process seemed to go to completion but I could not find the new, projected coverage anywhere. What is going on? I am baffled!


texts:
I'm trying to add popups for point icons using Leaflet.js and ran into this error: "Uncaught TypeError: points.eachLayer is not a function"

Here's my code:

function addPoints(data){
    for (var feature in data.features){
        if(data.features[feature].properties.activity == "hunting"){
            points = L.marker([data.features[feature].geometry.coordinates[1],data.features[feature].geometry.coordinates[0]], {
                icon: L.mapbox.marker.icon({
                    'marker-color': '#FF5A5C',
                    'marker-symbol': 'bus'
                })
            }).addTo(map).on({
                mouseover: function(e){
                    this.openPopup();
                }, mouseout: function(e){
                    this.closePopup();
                }
            });
        }
    };
    drawPopup(data);
} //end addPoints
function drawPopup(data){
    points.eachLayer(function(layer){
        var popupContent = "Park: " + data.features[feature].properties.name + "<br />" + "Activities: " + data.features[feature].properties.activity;
        layer.bindPopup(popupContent);
    });
}

I've used nearly identical code before and have never run into this problem, what am I doing wrong?


texts:
is it possible to integrate the following WPS-LINK in QGIS?
http://gdz.bkg.bund.de/wps_koordtrans?REQUEST=GetCapabilities&SERVICE=WPS


texts: Error message: Value cannot be null. Parameter name: symbolData

texts:
I am trying to reference a web map in Geocortex Essentials but it will not let me and gives the following error message: Value cannot be null. Parameter name: symbolData.

I isolated the layers in my web map and see that there are two layers causing this issue. The layers that don't work and the ones that do show no schema differences (they have different features). But no feature has a null value and I tried changing the symbology of the problematic layers to match the ones that work. Still no luck. Any ideas?


texts:
I need to save my vectors layer back to PostGIS, which if I read correctly I can do using GeoServer, however I'm confused.  All the examples ask you to import a shapefile.  I don't have a shapefile.  This will be a empty vector layer on top of a couple of base layers which the user will add polygons, points etc. etc.

How do I create the tables to hold the vector data.  Also how do I know the structure of the tables.  e.g columns, column names, primary / foreign key, relationships etc. etc.


texts: Move legend if it overlaps features within dataframe using ArcPy

texts: Trying to find a way programmatically (arcpy) move the legend if it intercepts features within a data frame, in the scenario below, if the legend obscures the view of the AOI, then I want it to move to a different corner until its not an issue. This has to be on top of the data frame as opposed to making the data frame smaller and just putting it to the side.


texts: Dialog Table Name Prompt Error In MapBasic_Mapbasic

texts:
I am trying to prompt the Table name from the opened tables and to use it in the rest of the program, but what i get is (Table (i) is not open) after clicking ok
the Sub Procedure is below, please help

Sub Selection

Dim num_open_Tabs, Ix  As Integer
Dim Tab_names(50) As string
Dim Sector_table as string
Dim Sector_Table_Name as String
Dim sector_Table_Obj as Object

Win_ID = FrontWindow()
num_open_Tabs = NumTables() 'Create Array of the Layers Names

Redim Tab_Names(num_open_Tabs)
For Ix = 1 to num_open_Tabs
   Tab_names(ix) = TableINFO(Ix, Tab_INFO_NAME)
Next

Dialog 'Prompt the wanted Layer name
   Title "Choose Cell_Data Table"
   Control StaticText
      Position 10,10 Width 250 Height 10
      ID 1
      Title "Please Choose The Table Containing Your Sector Info"

   Control popupmenu
      Title From Variable Tab_names
      Position 10,20 Width 250 Height 15
      Id 2
      Value Sector_table
      Into Sector_table
   Control StaticText
      Position 10,50
      Title "Enter Buffer Width"

   Control EditText
      Position 10, 60
      Value "200"
      Into Sel_Width
   Control OKButton
   Control CancelButton

If CommandInfo (CMD_INFO_DLG_OK) then
   Sector_Table_Name = TableInfo(Sector_Table, TAB_INFO_NAME)
Else
   Exit Sub
End If


texts: Installing ArcGIS 10.1 on SQL Server 2008 R2 Cluster

texts:
We're running into some performance issues when accessing large enterprise geodatabases on a SQL Server 2008 R2 cluster.  The database itself performs well and returns queries fine but presenting the information back to ArcGIS for Desktop 10.1 is prohibitively slow.

Connection is direct connect.
Data type in SDE is DBO.
Machines are VMs with 16Gb RAM.

Are there any tips or tricks or known documents for installs of these on SQL Server clusters where a node instance is presented?


texts: How to calculate the optimum bus route using pgrouting?

texts:
My setup for my web page is OpenLayers, Postgresql/postgis and pgrouting.

What I do not understand is how I am going to do the routing.I have read many tutorials but I don't understand what my next step will be.

So I have downloaded an shp file of osm data into my computer and imported the 'roads.shp' into my database. In the database, I have also all of my bus route stops as tables of points.

How do I use pgrouting from now on? What algorithm to use and how to create the topology? Do I need to merge the roads shp with the tables of the bus route stops ?

For additional info please don't hesitate to ask!


texts: GeoServer 2.9.2 does not update Bounding Box after layer shapefile is generated second time

texts:
ImageMosaic store was created for raster layer. Everything OK.

Now additional .sid files were added so new shape was generated manually with gdaltindex command. When new shapefile is created I go to layer configuration in GeoServer and press "Compute from data" but BBox does not change (I know for sure it now has larger BBox). If I do a GeoServer restart and then press "Compute from data" for this layer only then BBox gets recalculated. That is an issue one.

Even when I successfully restart GeoServer and recalculate BBox issue two persists: when opening GWC GUI for layer caching I see that BBox still has the old values. To fix this I have to go to ../data/gwc-layers/LayerInfoImpl-xxxx and manually change coordinates.


texts: Using ArcGIS to assign values from a shapefile to a new field of a second shapefile if they are intersecting

texts:
I have a polygon feature set (blue) which has the shape of a raster dataset and I have another shapefile (green) where the polygons are seperated by the "borders" of the raster polygon file. What would be the fastest way to assign the values of one column of the green shapefile to a new column in the blue shapefile if they intersect? Also if they do not intersect assigning a value like 9999 would be perfect. So far I got the best result with union but I would like avoid creating a new file and I also don't need every column from the green shapefile.

My files look like this:


texts: PLSS shapefiles - alternative to geocommunicator.gov?

texts: The BLM's website geocommunicator.gov is down and I'm wondering if anyone knows of a second source for PLSS shapefiles.


texts:
How can I make marker Icon appear for each point. I added the geoJSON with 65 points and have all been validated. However the icons appeared as a broken png.
I proceeded to add a png
I called in an icon from the img folder : "img/marker.png"
I passed the icon function L.geoJson:  //this is where the icon does not appear for every single point but only one big icon appeared.

Can you guys check the code for why the marker does not appear for all the points but for just one

Find here all the files.
https://github.com/saviour123/malariaRes

//JS and Leaflet Stuff
//creating the map with div
var map = L.map('map', {
    center: [5,28],
    zoom: 3,
    minZoom: 2,
    maxZoom: 18
});

//add tiles
L.tileLayer('http://{s}.tiles.mapbox.com/v3/ianmule.bg2v5cdi/{z}/{x}/{y}.png',{attribution: "Mapbox"}).addTo(map);

//create map icon
var marker = new L.Icon({
    iconUrl: 'img/marker.png',
    iconSize: [38, 95],
    iconAnchor: [22, 94],
    popupAnchor: [-3, -76]
});

//the function that binds the icon to the layer
function myIcon(feature, layer) {
    layer.setIcon(marker);
};

L.geoJson(fielddata,{
    onEachFeature: myIcon
}).addTo(map);


texts: Lat-Long values get rounded off when stored in a column type Geography in postgres

texts:
I have some lat/long coordinates eg:

  (64.123456789123456, -110.12341234123412)

I store them in the table using a query that looks like this:

  INSERT INTO myTable (name, geo) VALUES ('test', ST_GeogFromText('SRID=4326;POINT(-110.12341234123412 64.123456789123456)') );

but when I try to retrieve the lat-long back from the table, they seem to be rounded off, in this case, I get

  lat: 64.1234567891234

and

  long:-110.123412341234

  select ST_X(geo::geometry) as long, ST_Y(geo::geometry) as lat from myTable;

I am probably not using the right data types resulting in loss of values, but haven't been able to figure out.

I am new to PostgreSQL.

How do I not lose the coordinate precision?


texts: Select by Expression in QGIS Graphical Modeler

texts:
I am trying to use select by expression or extract by expression in QGIS Graphical Modeler.
I add a vector layer for input (have also tried vector features) I add the algorithm, select the input layer, then go to build the expression and there are no fields listed under Fields and Values.
If I manually input the expression I get an error that the column is not found.


texts:
I am trying to have multiple labels on a non-geographic map (image), but when I zoom in and out, I would like the text to grow and shrink in size corresponding to the zooming.

Here is what I have:

JS

 var mapExtent = [0.00000000, -16384.00000000, 16384.00000000, 0.00000000];
 var mapMinZoom = 0;
 var mapMaxZoom = 6;
 var mapMaxResolution = 1.00000000;
 var mapMinResolution = Math.pow(2, mapMaxZoom) * mapMaxResolution;;
 var tileExtent = [0.00000000, -16384.00000000, 16384.00000000, 0.00000000];
 var crs = L.CRS.Simple;
 crs.transformation = new L.Transformation(1, -tileExtent[0], -1, tileExtent[3]);
 crs.scale = function(zoom) {
   return Math.pow(2, zoom) / mapMinResolution;
 };
 crs.zoom = function(scale) {
   return Math.log(scale * mapMinResolution) / Math.LN2;
 };
 var layer;
 var map = new L.Map('map', {
   maxZoom: mapMaxZoom,
   minZoom: mapMinZoom,
   crs: crs
 });

 layer = L.tileLayer('{z}/{x}/{y}.png', {
   minZoom: mapMinZoom,
   maxZoom: mapMaxZoom,
   noWrap: true,
   tms: false,

 }).addTo(map);
 map.fitBounds([
   crs.unproject(L.point(mapExtent[2], mapExtent[3])),
   crs.unproject(L.point(mapExtent[0], mapExtent[1]))
 ]);
 L.control.mousePosition().addTo(map)

 var yx = L.latLng;
 var xy = function(x, y) {
   if (L.Util.isArray(x)) { // When doing xy([x, y]);
     return yx(x[1], x[0]);
   }
   return yx(y, x); // When doing xy(x, y);
 };

 L.marker(xy(-5213, 57.0), {
   icon: L.divIcon({
     className: 'text-labels',
     html: 'Testing a long sentence here'
   }),
 }).addTo(map);
 L.marker(xy(-3213, 57.0), {
   icon: L.divIcon({
     className: 'text-labels',
     html: 'Image 2'
   }),
 }).addTo(map);
 L.marker(xy(-5213, -1057.0), {
   icon: L.divIcon({
     className: 'text-labels',
     html: 'Image 3'
   }),
 }).addTo(map);
 L.marker(xy(-10213, 3057.0), {
   icon: L.divIcon({
     className: 'text-labels',
     html: 'Image 4'
   }),
 }).addTo(map);
 L.marker(xy(10213, -5057.0), {
   icon: L.divIcon({
     className: 'text-labels',
     html: 'Image 5'
   }),
 }).addTo(map);

CSS

html,
body {
  width: 100%;
  height: 100%;
  margin: 0;
  padding: 0;
}

#slider {
  position: absolute;
  top: 10px;
  right: 10px;
  z-index: 5;
}

#map {
  position: absolute;
  top: 0;
  right: 0;
  left: 0;
  bottom: 0;
  top: 0;
  z-index: 0;
  background: none;
}

.text-labels {
  font-size: 1em;
  font-weight: 700;
  cursor: grab;
  white-space: nowrap;
}

Here is a fiddle


texts: Rescaling and reoffsetting a point cloud with lidR?

texts:
I'm trying to analyse a a 3 x 3 metres scan derived from terrestrial photogrammetry (i.e. local coordinates, with 0,0 in the centre of the point cloud) to obtain a DSM.

Importing the file in lidR, I get the following warnings:

Warning messages:
1: Invalid header: X scale factors should be factor ten of 0.1 or 0.5 or 0.25 not 2.23339664936066e-09
2: Invalid header: Y scale factors should be factor ten of 0.1 or 0.5 or 0.25 not 2.33016455173492e-09
3: Invalid header: Z scale factors should be factor ten of 0.1 or 0.5 or 0.25 not 4.71837997436523e-10

I specified the scale to 0.001 using lasrescale (suggested here too).

Then, I tried using lasreoffset to shift the point cloud to the correct geographical coordinates (here a similar implementation using LAStools), but I get this error:

las = lasreoffset(las, xoffset = 1306977.096, yoffset = 5661354.023, zoffset=265.8018)
#> Error: Incorrect xoffset: integer overflow.

When manually adding the offset to each point (e.g. las@data$X + 1306977.096, etc), the procedure works.

Finally, once I try to rasterize the point cloud using:

r = grid_canopy(las, res=0.01, dsmtin())

I get another error:

Error: Internal error in C_interpolate_delaunay: xy coordinates were not converted to integer. Scale factors are likely to be invalid.

It looks like I'm using the wrong tools for the purpose, or totally missing the point. Any suggestions?

Furthermore, can anybody help me understanding what does "integer overflow" refer to in this case?


texts: Symbolising with one expression and labelling with another using QGIS?

texts:
Edit - I forgot that I had posted this until I was searching for information on this.

I have 2 land type layers (soil, vegetation types etc), one of the layers shows the land types for entire paddocks while the second layer shows only the land types in the paddock within a certain distance of water. What I want to do is have both layers with the same symbology but displaying different labels containing the land unit, land unit description and the total area of the land unit. Obviously the total area would be different between layers.

Further edit:

As an example based on the picture below, I would like to catagorise based on the first part of the text string eg "De" or "Bs-1" so that all layers have the same symbology regardless of their different areas.

I would then like to create the category descriptions with an expression similar to: 'MAPUNIT' || " - " || 'DESCRIPTION' || round(sum('AREA', 'MAPUNIT'), 2))

I'm wondering if it is possible to categorise layers in QGIS using one expression but label it using another? I have several layers based on the same dataset but clipped so they have different areas. I want to categorise them so the symbology matches between the layers but label the categories with the areas as well.

At the moment I am using the same expression for each layer but have to manually change the category colours to match.

In ArcMap, I used the advanced symbol legend label tool from maplogic for the this.

An example of what I'm doing. The colours and description stay the same but the area change based on buffered areas. I'm using an expression to categorise the layers e.g. 'MAPUNIT' || " - " || 'DESCRIPTION' || round(sum('AREA', 'MAPUNIT'), 2))


texts:
Solution:

The limitations of the Esri + Oracle (per 12c) were the cause. 12c Oracle will allow for Multi-Tenant DB, allowing us to have multiple SDE GeoDB as we need. This will meeting the need for one owner, and multiple creators in a schema.

Background:  Oracle 11g, ArcServer 10.3.1

Our group is trying to organize our geodatabase based on application, much like how SQL Server creates its instances.

All our tables are at the master SDE level, which means that there are impacts on other application during schema deployments or esri upgrades.

We looked at User Geodatabase (user schemas), but they only allow for only editor/creator at the Geodatabase level.

Outside of creating one Oracle instance per application, what are others doing in order to allow for multiple owners and editors?


texts: Convert AutoCAd file (.dwg) to shapefile

texts:
I have a CAD file to convert into shapefile. Unfortunately, I don't have access to AutoCAD, MapInfo or ArcGIS atm. All I have access to is QGIS.

I'll need to convert this CAD file to shapefile to be used in a specialised software which can use/edit shapefile.

I tried couple of online tools without any success. Any help? BTW the size of the CAD file is less than 2MB.


texts: How do I precisely select lines within a polygon in QGIS?

texts:
I want to select lines (bahnstrecken) within a polygon (bundesland). My goal is a selection shown in the following screenshot, which I did manually. Before I used two tools: I converted the polygons to lines and splitted the linies (bahnstrecken) with the new line-layer which was formerly a polygon. The result is a lines layer, which geoemtries end exactly at the polygon borders.

The problem is this: When I use different spatial query methods (intersects, contains, touches, fully contains) I never get the result I want. How can I reach the selection result shown in the first screenshot above? I think I should use the WITHIN method, but the results of the spatial selection tool appear quite strange to me:
INTERSECT:

TOUCHES:

OVERLAPS:

WITHIN:


texts: How to handle with "no-Geometry" shapes in FME?

texts: I have many shape files that I am reading on FME (points, lines and polygons), editing them little bit and writing back the same format. However I am having an issue that I as far as process comes to an empty shape file while reading, process stucks and doesn´t continue enymore (though no error message pops up)! Is there a way that can ignore those empty shape files without any geometry (without any feature in it)?


texts: is there significance difference between qgis version 1.8.0 & 2 to install plugin?

texts: I am a Qgis learner on widows xp and I want to know whether I can install semiautomatic classification plugin onto Qgis 1.8.0 in the same way for Qgis 2.  is there significance difference between  the two version in accepting the plugin?


texts: What is the use of Automatic overview support in Qgis 2.0?

texts:
Automatic overview support in map frame is one one of the new feature added in QGIS 2.0. In QGIS 2.0 Documentation explanation about this is given as:

Automatic overview support in map frame:

"Need to show the current area of the main map frame in a smaller overview window. Now you can. The map frame now contains the ability to show the extents of other and will update when moved. Using this with the atlas generation feature now core in the composer allows for some slick map generation. Overview frame style uses the same styling as a normal map polygon object so your creativity is never restricted".

I can't understand its usage clearly. Can anyone explain it briefly?


texts: Dataframe data handling or SQL Query in ArcMap

texts:
I'm trying to perform an SQL query in ArcMap with the following data set, and I have not been able to solve the problem.

The idea is to keep all the information of the columns, but for example if the column "sens" has the value "S1" discard the value "S2_adjusted" that has a value of "NaN" but leave the column "S1_Adjusted". The same process for the data with value "S2" in the "sens" column, that is, leave only the value of "S2_adjusted" and discard "S1_adjusted". Without losing information in other columns.

I am using Arcmap 10.6.1 and the data comes from an Excel file.

Do you think that an SQL query is the correct approach?


texts: How can I split a polygon by network line in QGIS?

texts:
I have a road network as line layer and a land cover as polygon layer, I want to find number of subpolygon for each polygon that divided by line layer in QGIS?
I used union function in Geoprocessing Tools in Vector menu of QGIS and then polygonize function to convert union result to polygon but in table of new polygon there is NULL data.


texts:
I am trying to read a Cloud Optimised GeoTiff (COG) from AWS S3. I can successfully run

gdalinfo /vsis3/cog-testing/HP_COG.tif

And see the details of the file and it says Layout: COG

But if I attempt to load the same file into QGIS 3.4 using the Data Source Manager, using the AWS S3 Protocol.

But when i try I get an error

Invalid Layer: GDAL provider Cannot open GDAL dataset /vsis3/cog-testing/HP_COG.tif: `/vsis3/cog-testing/HP_COG.tif' does not exist in the file system, and is not recognized as a supported dataset name. Raster layer Provider is not valid (provider: gdal, URI: /vsis3/cog-testing/HP_COG.tif

I am wondering if it is a region problem as my files are not the in the default region.

I have also tried using the standard HTTP protocol and also get an error.

Invalid Layer: GDAL provider Cannot open GDAL dataset /vsicurl/https://cog-testing.s3-eu-west-1.amazonaws.com/HP_COG.tif: `/vsicurl/https://cog-testing.s3-eu-west-1.amazonaws.com/HP_COG.tif' does not exist in the file system, and is not recognized as a supported dataset name. Raster layer Provider is not valid (provider: gdal, URI: /vsicurl/https://cog-testing.s3-eu-west-1.amazonaws.com/HP_COG.tif

And when I try the viscurl from my command prompt I get a 403 error.

gdalinfo /vsicurl/https://cog-testing.s3-eu-west-1.amazonaws.com/HP_COG.tif

Can anyone suggest what is going wrong?


texts: Adding a feature to Openlayers Vector source

texts:
I want to;

get the features from a vector source and
for some features, use their coordinates to construct new features, e.g. ol.geom.Circle buffers.

My sample vector source (epsg:3857 projection), with two features (points) that I would like to read their coordinates from:

var sampleSource = new ol.source.Vector({
    url: (basePath + '/data/scanner_locations.json'),
    format: new ol.format.GeoJSON()
    "EPSG:3857"})).readFeatures(basePath + '/data/scanner_locations.json')
});

This layer renders two points on my map.

For 1, to retrieve the source features, if I understand correctly you need to monitor the change event like this:

sampleSource.on('change', function(evt){
var source = evt.target;
if (source.getState() === 'ready') {

If it is 'ready' you have asynchronously retrieved the source, correct?

Following I try to add something to the source like so:

sampleSource.on('change', function(evt){
var source = evt.target;
if (source.getState() === 'ready') {
    var numFeatures = source.getFeatures().length; // if 0, quit ?
    var features = source.getFeatures();
    console.log(features); // gives me the two features
    features.forEach((feature) => {
        if (feature.values_.id === 1){ // if the id property of the feature is 1
            // feature.getGeometry().flatCoordinates should be the center of the circle
            source.addFeature(new ol.Feature({
            // coordinates for the circle center are hardcoded here, radius is 500 meter
                geometry: new ol.geom.Circle([486524.0550283753, 6801035.6531745335], 500),
                name: 'the buffer'
                // id: 500
            }));
        } // end if
    }) // end forEach
}});

This creates a huge stack size, so something is clearly going wrong.


texts: Set 00000 Input Transparent Color with GeoServer Rest API

texts:
I need to change parameters (Input Transparent Color) from layers in GeoServer. I have a lot of images to do it. So, I need to do this edit using GeoServer REST API.

I did the code below but it doesn't work.

curl -u admin:geoserver -v -XPUT -H "Content-type: text/xml" -d "<coverage><parameters><entry><string>*InputTransparentColor*</string><string>000000</string></entry></parameters><enabled>true</enabled></coverage>" http://localhost:$81/geoserver/rest/workspaces/sentinel2/coveragestores/$covareganame/coverages/$covareganame.xml

Is it possible to change layer's parameters from the curl commands?
What's wrong in my script it not change input transparent color to '00000'?


texts: Display time param as decoration in WMS request

texts:
Is there a way to add dynamic text decoration via a geoserver layout.xml that reflects the &time param of a geoserver WMS request?

I know from here how to display frame parameters as decorations for animations by using
value="${avalue}" in layout.xml and each animation frame shows the avalue param correctly, but for a non animation request there is no &avalue just &time and
value="${time}" doesn't seem to work.

Here is the layout that works for animations:

<layout>
  <decoration type="text" affinity="bottom,left" offset="16,6">
    <option name="message" value="Leaf Out Date: ${avalue}"/>
    <option name="font-size" value="22"/>
    <option name="font-family" value="Arial"/>
    <option name="halo-radius" value="2"/>
  </decoration>
  <decoration type="legend" affinity="bottom,right" offset="6,-6" size="0,0"/>
</layout>

My non animation request looks like this:

  http://myserver/geoserver/si-x/wms?service=WMS&request=GetMap&layers=si-x:average_leaf_ncep_states&format_options=layout:leafout_web&styles=si-x:leafout_bimonthly&time=2016-2-18&bbox=-125.020833333333,24.0625,-66.479166666662,49.937500000002&width=1400&height=600&srs=EPSG:4269&format=image/png


texts:
I am trying to detect green roofs in a city and would like to compute the total area of green roofs (also on a per roof basis).
I used a shapefile containing all roofs and used that to extract the roofs from an orthophoto and calculated NDVI through the raster calculator. Now I would like to compute the area of the green roofs. I'm hoping to get some kind of attribute table that shows me if a given building has a green roof and much m².
I tried using reclassify by table to create two classes (no vegetation & vegetation) and polygonized the output which sometimes resulted in hundreds of polygons within the boundaries of a single building.
I'm unsure about how to proceed or if I am even on the right track.
Please keep in mind that I am using QGIS.


texts: Enlisting all OpenStreetMap (OSM) street names for General Transit Feed Specification (GTFS) routes?

texts:
I have exported OpenStreetMap (OSM) data of a city to my local postgresql db using osm2pgsql.

I would like to add the General Transit Feed Specification (GTFS) routes information (say bus routes) into planet_osm_roads table as tags - say route_id(s) and agency_id(s).

All the open source GTFS tools seem to use gtfs stops as nodes. Opentripplanner, for example, does use osm street data while building the graph nodes to better understand the stops but the osm part stops right there.

I could probably iterate over all the routes, enlist all the stops in each of the trips and use routing library like PgRouting to determine the path (i.e street names) between those stops.

This seems to be a decent way, but there is no guarantee that the paths are indeed that used by the bus routes.

Any suggestions how should I be approaching this problem?


texts:
I am trying to delineate a watershed using r.water.outlet.

I have created a flow direction raster using r.watershed (threshold = 50) (cell size is 300m).

I use this flow direction raster in r.water.outlet
I then specify EXACT coordinates of a location on a stream network, also created using r.watershed (in wgs84 - crs that the grass porject is in).

a blank 'black' raster is returned. The style tab in properties has max:1. min:1 (when compute min max).

I have re-run the analysis several times with different output from r.watershed (i.e. using different thresholds) - but has the same result.


texts: GDAL - How to get geometry type of a MapInfo TAB file?

texts:
I am using the OGR MapInfo File driver to read a MapInfo TAB file. However, the function GetGeomType() returns 0 which means 'wkbUnknown` even if the features in the TAB are only multi-polygons.

Is it because MapInfo TAB can store many different GeomTypes in the same file?

If I take for granted that only one geometry type would always be in the TAB, how could I retrieve it?

        driver = ogr.GetDriverByName("MapInfo File")
        datasource = driver.Open(os.path.join(dirname,shapefileName))
        layer = datasource.GetLayer(0)
        geometryType = layer.GetGeomType()

        >>>geometryType
        >>>0


texts: Calculating how many of certain feature lie within various buffer rings?

texts:
I am using ArcMap and I have a series of real estate parcels mapped along with a series of wind turbines. I have created multiple buffer rings around each real estate parcel with distances (0.5, 1, 1.5, and 2) miles.

What I am trying to do is determine how many turbines fall within each buffer ring for EVERY real estate parcel. The near tool is almost what I need, but not quite. I am trying to get an output of something like this.

And so on.

Is there a simple tool that can make this happen?


texts: Communication with between webbrowser and QGIS desktop

texts:
I would like to send requests from a webbrowser to QGIS-desktop,
For example sending a WMS or WFS URL and add that in QGIS.

What options does it exist for this ?
I have in earlier systems doen this intgration through a database table that
the browser saves information into and that ArcGIS-desktop reads from.
That are though a quite Heavy solution.

Would it be possible to run QGIS-desktop as a tiny webbserver on a port like 8090 and let QGIS accept incling requests from the browser like
http://localhost:8090/addLayer?URL=http://www.mywms.se/aWMSServer

Using something like whats discussed in:

https://ruslanspivak.com/lsbaws-part1/

I have also previsouly tested in a .net application to start the webbrowser within a .net application as host. Then I could call back to parent application with windows.external in Javascript. But maybe that is only working for Internet Explorer and .Net desktop applications.

Anyone have an idea or better (modern) idea for communication between browser and a desktop app?

/Michael


texts: Unable to view cached map service unless logged in as admin in GIS Server

texts:
I recently re-created a cache for a map service due to SDE server changes. When the cache finished I was unable to view it unless I manually went to the cache folder and viewed each tile one by one in photoshop.

After a little research I discovered the virtual directory in the cache directory was incorrect. After changing the virtual directory I can now view the cached service in ArcCatalog through my admin login, but cannot view it through the "Use GIS Services". This happens on both the internal and external sites I have.

I've checked security/permissions on everything I can think of, folders, groups, etc. I'm almost certain this is an easy fix, but I'm coming up with nothing. Any help would relieve some stress!


texts: Python Script to convert multiple KMZs to single GDB

texts:
I am trying to extract all KMZ and kml files from a folder but it's giving me an error. File does exist in that folder.

    for kmz in arcpy.ListFiles('*.km*'):
        print ("CONVERTING: " + os.path.join(arcpy.env.workspace,kmz))
        arcpy.KMLToLayer_conversion(kmz, outLocation, "", "NO_GROUNDOVERLAY")

ERROR 000732: Input KML File: Dataset NMR102020202_XZCD7S_4850NW_ABCD5W_Poin.kmz does not exist or is not supported

I am not sure what's missing in this code.


texts: How to export a dataset with "SpatialPolygonsDataFrame" as a shapefile

texts:
How to export spatial datasets in R, such as tornados dataset in GISTools package, as a shapefile with the exact same spatial reference and the exact values?

library(GISTools)

data(tornados)

Desired:

"SpatialPolygonsDataFrame" to "shapefile"
"(Forced) SpatialPointsDataFrame" to "shapefile"
"(Non-Forced) SpatialPointsDataFrame" to "shapefile"
"SpatialLinesDataFrame" to "shapefile".


texts: Modify logo + copyright on Bing Maps REST API

texts:
I'm attempting to use Bing Maps REST API to track traffic flow over time.
To do so, I mosaic many 640x640 pixel tiles together, using this template
The issue is that when I mosaic all these tiles, the Bing map logo and API obscure many parts of the map, like this:

I know that it is possible to modify the logo with Bing Maps V8 Web Control by accessing the MapOptions Object , but I'm not sure how to access these with the API Rest (newbie with REST services over here).

For the record, I am aware of the terms and conditions of use for the Bing REST API, and am not planning to display these maps without a logo (I just would rather not have 50 logos/copyright info per map).