texts: Displaying flow direction using arrow in ArcGIS Desktop?

texts:
I created a flow direction raster out of a DEM using ArcGIS Spatial Analyst Tools/ Hydrology, so the directions are coded in raster values as displayed below.

[http://desktop.arcgis.com/de/arcmap/10.3/tools/spatial-analyst-toolbox/flow-direction.htm]

Then, changing the layer style in the Layer Properties to an arrow, I cannot set the right orientation as coded by the ArcGIS flow direction tool (first image). So as you see in the second image, the arrows don't indicate the right direction, which should actually lead to the sink (blue polygon).

Instead of using the layer properties tried the "Utility Network Analyst Toolbar", but the flow direction raster layer cannot be selected.

Can anyone help me where to adjust the orientation of the arrows?


texts:
How can I get current coordinate system type from opened file?
I tried
 var coord = ESRI.ArcGIS.Geometry.esriSRProjCSType.... ???;
Can you tell me, what should I do?

UPD:
I need to save current map with all layers and geometry into a shapefile. That's why I have to get current coordinate system.


texts: Creating polygon feature class from area in one not covered by another using ArcGIS Desktop?

texts:
I am working on a project doing spatial analysis at a statewide level. I am looking at a form of transportation coverage in a state, and part of my analysis involves buffers and drive time polygons. I am using a geodatabase (therefore ArcGIS) for all my work, so new features will be exported as feature classes in feature datasets.
To visualize: I have said state border, and some dissolved buffer rings around some points, and they cover parts of the state. I want to create a feature class for the parts of the state these buffer rings don't cover, AKA gap areas.
Is there a tool I can use or create to do this?
I know I could draw them in manually in editor, but I am unsure if that's the most precise way to generate these feature classes. To me, I am thinking this operation would be like an inverted "Clip" in that I take the two layers, and instead of pulling out the clipped area, it would take everything not in the clipped area.
I am using a Desktop Basic license.


texts:
I have a cached map service that was published on ArcGIS Server 10.3. The tiles were built with no error. Using the ArcGIS/Bing Maps scale, I set the tiles to build to 1:577,791, I also set the extent of the layer to turn off in beyond 1:577,791. I zoom into the 1:577,791 scale, the layer is still checked on but no tiled data draws. When I zoom into the next scale. 1:288,895 the cached layer turns off/is grayed out, and no data appears, as it should.

Any thoughts as to why there is no data drawing at the 1:577,791 scale? I would like the tiles to draw at that scale.


texts: histogram of classified image in Google Earth Engine

texts:
I'm using Google Earth Engine to run a supervised classification and I want to visualize, for each class in my classification, the distribution of DN's or Reflectance in each band of a Landsat image.

The way i visualize the result is as follows:
I would like to have a histogram like the one below for each class in my classified image

How do i select a property from a layer (in this case cand cover type to create a histogram of those pixels?

using GEE's examples:

    // Upsample MODIS landcover classification (250m) to Landsat
// resolution (30m) using a supervised classifier.

var geometry = /* color: #ffc82d */ee.FeatureCollection(
        [ee.Feature(
            ee.Geometry.Polygon(
                [[[29.9761962890625, 31.592573630393357],
                  [29.981689453125, 30.111869849235248],
                  [32.574462890625, 30.116621582819374],
                  [32.530517578125, 31.620643692450585]]]),
            {
              "system:index": "0"
            })]);

// Use the MCD12 land-cover as training data.
var modisLandcover = ee.Image('MODIS/051/MCD12Q1/2001_01_01')
  .select('Land_Cover_Type_1');

// A pallete to use for visualizing landcover images.
var landcoverPalette = [
    'aec3d4', // water
    '152106', '225129', '369b47', '30eb5b', '387242', // forest
    '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40',  // shrub, grass, savanah
    '111149', // wetlands
    '8dc33b', // croplands
    'cc0013', // urban
    '6ca80d', // crop mosaic
    'd7cdcc', // snow and ice
    'f7e084', // barren
    '6f6f6f'  // tundra
];

// A set of visualization parameters using the landcover palette.
var landcoverVisualization = {palette: landcoverPalette, min: 0, max: 17, format: 'png'};
// Center over our region of interest.
Map.centerObject(geometry, 11);
// Draw the MODIS landcover image.
Map.addLayer(modisLandcover, landcoverVisualization, 'MODIS landcover');

// Load and filter Landsat data.
var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1')
    .filterBounds(geometry)
    .filterDate('2000-01-01', '2001-01-01');

// Draw the Landsat composite, visualizing true color bands.
var landsatComposite = ee.Algorithms.Landsat.simpleComposite({
  collection: l7,
  asFloat: true
});
Map.addLayer(landsatComposite, {min: 0, max: 0.3, bands: ['B3','B2','B1']}, 'Landsat composite');

// Make a training dataset by sampling the stacked images.
var training = modisLandcover.addBands(landsatComposite).sample({
  region: geometry,
  scale: 30,
  numPixels: 1000
});

// Train a classifier using the training data.
var classifier = ee.Classifier.cart().train({
  features: training,
  classProperty: 'Land_Cover_Type_1',
});

// Apply the classifier to the original composite.
var upsampled = landsatComposite.classify(classifier);

// Draw the upsampled landcover image.
Map.addLayer(upsampled, landcoverVisualization, 'Upsampled landcover');

// Show the training area.
Map.addLayer(ee.Image().paint(geometry, 1, 2), null, 'Training region');


texts: Trying to use layer saved as .lyr file with arcpy.mapping.Layer function results in error?

texts:
I am fairly new to the ArcGIS development world. I am attempting to write a script that will have an end result of selecting features by proximity to a selected portion of a different feature.
So far, my code executes the following steps (with varying degrees of success):

Create a layer from selected features,
Save the layer onto a disk location,
Add that layer back to an ArcMap document,
And then use the layer in a Geoprocessing tool (Select Layer By Location).

The reason I am attempting to do it in this order is that I am running my code from a python window external to the python window within an ArcMap document (i.e. the ArcMap document could be closed while I am running this code).
Thus far, I have been able to create a temporary layer from a selection of an original layer, and save a copy of this temporary layer to a file geodatabase (I used the arcpy.SaveToLayerFile_management function to save as a .lyr file). When I look at the file geodatabase through the Windows file explorer, the temporary layer is successfully saved there. However, when I look in the ArcMap catalog, the .lyr file is not visible. On top of that, any time I try to use the arcpy.mapping.Layer function to call the layer, I get the error:

Traceback (most recent call last):
  File "H:/Work/Python_Project/ExportGBSSDataV1.0/Spatial_Location_Test.py", line 46, in <module>
    addLayer = arcpy.mapping.Layer(database + '\\' + temp_CSlyr + '2.lyr')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\arcobjects\mixins.py", line 429, in __init__
    super(LayerMixin, self).__init__(lyrfile)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\arcobjects\_base.py", line 47, in __init__
    for arg in args))
ValueError: Object: CreateObject Layer invalid data source

What is going on here, and how do I fix it?

I’ve read other posts where the same error is encountered, but they were slightly different problems.

My Script is below:

# declaring initial variables and workspaces
arcpy.env.workspace = 'G:\\foo\\GDB.gdb'

database = 'G:\\foo\\GDB.gdb'

LRS_Loc = 'G:\\foo\\OtherGDB.gdb' #this layer was from a different Geodatabase

# Input for attribute to select
CS = raw_input('Enter the Highway Number and Control section with no spaces: ')

# create whereclause for MakeFeatureLayer_management function
CS_whereclause = '"ROADNAME" LIKE ' + "'" + str(CS) + "%'"

# Create variables for Layer names
temp_CSlyr = 'temp_CS'
BHLOClyr = 'Boreholes'

# Open mxd to work with Geoprocessing tools later
mxd = arcpy.mapping.MapDocument('G:\\bar\\Map.mxd')

#Create temporary layer from selected features
arcpy.MakeFeatureLayer_management(in_features=LRS_Loc + "\\HIGHWAYS",
                                  out_layer=temp_CSlyr,
                                  where_clause=CS_whereclause,
                                  workspace=database)

# Save temporary layer to file Geodatabase
arcpy.SaveToLayerFile_management(temp_CSlyr, temp_CSlyr, "ABSOLUTE")

# Add Layer to map (code stops working here)
df = arcpy.mapping.ListDataFrames(mxd)
addLayer = arcpy.mapping.Layer(database + '\\' + temp_CSlyr + '2.lyr')
arcpy.mapping.AddLayer(data_frame=df,
                       add_layer=addLayer)

# Select Layer By Location Geoprocessing tool
arcpy.SelectLayerByLocation_management(in_layer=BHLOClyr,
                                        overlap_type='WITHIN_A_DISTANCE',
                                        select_features=temp_CSlyr,
                                        search_distance='50' + " Meters")
del mxd, addLayer


texts: How to calculate polygon areas in acres or square miles instead of feet?

texts:
I have a layer in State Plane NAD83. When I use "Vector|Geometry Tools|Add/export geometry columns" or Field Calculator they both create the area in feet. What if I want acres or square miles? Thanks

How to calculate polygon areas in QGIS?


texts: How to merge groups of touching lines to multilinestring

texts:
I would like to join groups of lines into multilinestrings if they touch. In the plot below, I should get two objects of sf class MULTILINESTRING from the 6 lines. What is an easy way of doing this (when having 2 mio lines)? Is there maybe a pre-defined command that let's one do this?

data:

library(sf)
lines <- st_sfc(list(st_linestring(rbind(c(0,0), c(0.5, 0.7))),
                 st_linestring(rbind(c(0.5, 0.7), c(1, 1))),
                 st_linestring(rbind(c(2,0), c(1,1))),
                 st_linestring(rbind(c(1,1), c(1, 3))),
                 st_linestring(rbind(c(3,2), c(4,2))),
                 st_linestring(rbind(c(3,2), c(5,3)))))

I tried the following

library(igraph)
components(graph.adjlist(st_touches(lines)))

graph.adjlist transforms the list of elements that touch each other to a graph object and the function components identifies the components of each cluster.

But is there a way to do connect the multiple LINESTRING into MULTILINESTRING if they touch?


texts: Zoom to country with geoserver and openlayers

texts:
I would like to preposition an OpenLayers map on a country.
This is, when my page loads, the map should be positioned on the country with an appropriate zoom level. This should work for a (large) list of countries.

I have a GeoServer with a shapefile layer containing country polygons. The polygons have ISO country code attributes. Another layer shows a satellite background image.
When generating the country page on which the map is shown, I only know the country name (or ISO code), but not its coordinates and polygon details. (GeoServer and webserver are different machines)

My knowledge of OpenLayers and GeoServer is limited but if I understand correctly OpenLayers requests a certain map position and zoom level from GeoServer. I would thus need to retrieve position and zoom level from GeoServer in a first step, and then request this part of the map in a second one. This is if there is not a way to tell OpenLayers to zoom in on a feature I can identify by its attribute.

Any code example would be very much appreciated, but conceptual clarification would be equally valuable.


texts: Leaflet popups - preserving user input on close/reopen

texts:
I am working on a Leaflet map that allows a user to open a popup by clicking on a geometry and then generate dynamic content within the popup.  However, I find that on closing and then reopening the popup, the dynamic content is lost - I can think of a very lengthy workaround to preserve the content, but I'd like to see if anyone has a suggestion about a more direct way to do this, perhaps using Leaflet's methods... is there any way to preserve a popup's dynamic content when it has been closed and reopened?

Here is a slightly edited version of my code...

map.on('draw:created', function(e){

    //SETS THE BASIC FORMAT, INCLUDING BUTTONS FOR ADDING C
    popUpFields = "<button class='add_content'>Add content</button>

<button class='popup_save'>Save</button>";

    $(document).on("click", ".add_content", function(){
        $(this).next(".content_container").append("Hello popup!");
    });

    //ADDS THE LAYER TO THE MAP
    featureGroup.addLayer(e.layer);
    //AUTOMATICALLY OPENS THE POPUP
    e.layer.bindPopup(popUpFields).openPopup();

    e.layer.on('click', function(){
        //IF THE POPUP'S CONTENT HAS BEEN 'SAVED' AS AN ATTRIBUTE OF THE GEOMETRY - IN OTHER WORDS, THE USER HAS EDITED THE POPUP CONTENT
        if(e.layer.hello){
            //THIS IS WHERE A REOPENED POPUP NEEDS TO APPEAR WITH ALL OF THE CONTENT THAT HAD BEEN 'SAVED' - CODE BELOW DOES NOT WORK (OPENS WITH ONLY popUpFields CONTENT, NO ADDED CONTENT, BUT IS ILLUSTRATIVE OF WHAT I'M LOOKING FOR
            e.layer.openPopup();
        }
        else{
            e.layer.bindPopup(popUpFields).openPopup();
        };
    });

    $('.popup_save').click(function(){
        //THE USER SAVES THEIR POPUP CONTENT, WRITING THE ADDED CONTENT TO AN ATTRIBUTE OF THE LAYER
        e.layer.hello = $(this).prev().html();
    });

});

So in summary, what happens here is that a user adds a geometry, a popup opens, the user clicks the 'Add content' button a few times, appending a few strings to a div in the popup... the user saves those strings as an attribute of the geometry and then closes the popup to do other stuff.  Later, they wish to reopen the popup to add more strings, but all of the original added strings have been lost from the view.  I would like it so that the added content is not lost when the popup closes.  I actually have implemented a bit of a 'hack' solution that rebuilds the popup with the dynamic content every time it is reopened, but I feel that there must be a more elegant solution, where the state of the popup is preserved between closings and reopenings rather than freshly rebuilding it every time.  Is there?


texts: qgis-web-client stuck on "Loading map" on some maps

texts:
Some maps work for me to load, while some of them never do.
They do load the layer names, but nothing more.
The javascript console in Chrome returns this:

Uncaught TypeError: Cannot read property 'width' of undefined WebgisInit.js:450postLoading WebgisInit.js:450Ext.tree.AsyncTreeNode.listeners.load WebgisInit.js:236EXTUTIL.Event.fire ext-all-debug-w-comments.js:518EXTUTIL.Observable.fireEvent ext-all-debug-w-comments.js:159Ext.data.Node.Ext.extend.fireEvent ext-all-debug-w-comments.js:54140Ext.extend.loadComplete ext-all-debug-w-comments.js:55372(anonymous function) ext-base-debug.js:969Ext.extend.processResponse QGISExtensions.js:256Ext.extend.handleResponse ext-all-debug-w-comments.js:56415Ext.extend.handleResponse ext-all-debug-w-comments.js:6528handleTransactionResponse ext-base-debug.js:1887checkResponse ext-base-debug.js:1935(anonymous function) ext-base-debug.js:922

And if the project file (.qgs) contains a special character like åä or ö, it does not even load the layernames, and outputs this:

Uncaught TypeError: Cannot read property 'request' of undefined QGISExtensions.js:196Ext.extend.processResponse QGISExtensions.js:196Ext.extend.handleResponse ext-all-debug-w-comments.js:56415Ext.extend.handleResponse ext-all-debug-w-comments.js:6528handleTransactionResponse ext-base-debug.js:1887checkResponse ext-base-debug.js:1935(anonymous function) ext-base-debug.js:922

It does not matter if the shp-files themselves contain special characters, as I tried that on a working map without any problems.


texts: ArcMap ModelBuilder: problem with %name% variable when used in iteration for projection

texts:
I have prepared a model which should project all shape files ( currently in MUTM 81)  in a folder to MUTM 84. But when I use “%Name%” in output dataset for geodatabase as “D:\Soil Test\2881\D15\Dhuduwa.gdb\%name%”  I get message

 “Executing (Iterate Feature Classes): IterateFeatureClasses "D:\Soil Test\2881\D15\test" # # RECURSIVE
Start Time: Mon Mar 02 16:58:01 2015
Succeeded at Mon Mar 02 16:58:01 2015 (Elapsed Time: 0.17 seconds)
Executing (Project): Project "D:\Soil Test\2881\D15\test\BUILD_PT Point.shp" "D:\Soil Test\2881\D15\Dhuduwa.gdb\BUILD_PT Point" PROJCS['MUTM84',GEOGCS['GCS_Nepal_Nagarkot',DATUM['D_Nepal_Nagarkot',SPHEROID['Everest_Adjustment_1937',6377276.345,300.8017]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',84.0],PARAMETER['Scale_Factor',0.9999],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]] # PROJCS['MUTM81',GEOGCS['GCS_Nepal_Nagarkot',DATUM['D_Nepal_Nagarkot',SPHEROID['Everest_Adjustment_1937',6377276.345,300.8017]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',81.0],PARAMETER['Scale_Factor',0.9999],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]
Start Time: Mon Mar 02 16:58:01 2015
ERROR 999999: Error executing function.
The table name is invalid.
The table was not found. [BUILD_PT Point]
Failed to execute (Project).
Failed at Mon Mar 02 16:58:01 2015 (Elapsed Time: 0.06 seconds)

 “

I have 8 .shp files in that folder containing points, polylines and polygons.  What should I do?


texts: Alternatives to Add Surface Information and Extract Values from points

texts:
I am trying to add elevation data to an attribute table and am unable to use Add Surface Information and Extract Values as I have neither the spatial analyst nor 3D analyst license.
Are there any other tools/methods than can help me to achieve my goal?


texts: From Landsat spectral radiance (W/m2*sr*µm) to Global horizontal irradiance (W/m2)?

texts:
We are working in a solar energy project to locate potential sites for solar panel farms. We are using the Landsat 7 panchromatic band to derive radiance with GRASS using i.landsat.toar (W/m2*sr*µm units). Now, we have access to some ground data from sensors measuring global horizontal irradiance (GHI) which output is in W/m2 units.  I would like to compare these datasets for validation.

Is there a way to convert radiance in W/m2*sr*µm to irradiance in W/m2?

I have read in [1] that radiance = irradiance / π, but I don't completely understand, also, this doesn't take the spectral part (µm) into account.

[1] http://earth.esa.int/landtraining07/D1LB3-Su.pdf


texts: Seeking ArcPy script for choosing objects from table at random?

texts:
I am trying to write an ArcPy script that will select X number of rows from a table with N number of objects in it, then output that selection to a new layer or shapefile.  I need to be able to specify the number of selections I need at the time I run the script.  I am using ArcGIS 10.2 and working with a point FC.  Does anyone have any advice?

#import arcpy and set workspace
import arcpy
arcpy.env.workspace = "G:/583/Final/CRL_WaterSurfaces.gdb"

# set feature class and fields variables
fc = "MnEmbPoints"
fields = ["OID@"]

# list all OIDs in data set
with arcpy.da.SearchCursor(fc, fields) as cursor:
    for row in cursor:
        print("{0}".format(row[0]))


texts:
I’m trying to use the vertex tool in QGIS3 on a Mac to edit a shape file someone else sent me. I can edit points without any problems, but when I try to edit lines or polygons the tool doesn’t work. When I move the mouse close to or into a polygon, the vertices and line segments don’t get highlighted as they are apparently supposed to.

Why is this happening?


texts: Calculate distance of a point to the nearest road out of selected roads

texts:
I have a set of points that represent animal-vehicle collisions and a a layer of roads.
I want to calculate the distance between each collision point and the nearest road, but need to exclude the road where the collision happened. I am aware that the Near tool can help me with this, but can't find a way to exclude the road where the collision took place.

Any suggestions?


texts: Issues when exporting/zonal statistics rasters from .tif and .nc files

texts:
I am trying to aggregate at a 0.1º cell level, using zonal statistics, two variables that come from two different raster files (one in .tif, other in .nc format, links here and here).

Both seem to work fine in the R raster environment as, by loading and plotting them with

library(ncdf4)
r.works <- raster::raster('MOD13Q1_2000-02-18.250m_16_days_NDVI.tif')
r.dwork <- brick("spei01.nc") #this is a multilayer stack

rasterVis::levelplot(r.works)
rasterVis::levelplot(r.dwork[[1]]) #plot just the first layer

I obtain clear graphs, respectively attached. However, if I use the ZonalPipe function to do fast zonal statistics within a grid, I obtain a lot of NA values for the latter (not shown/plotted here for neatness), even though both rasters cover 100% of the grid of interest. I obviously set both rasters' crs equal to the grid's one before doing zonal stats. I thus investigated both objects and they seem to differ

> r.works
class       : RasterLayer
dimensions  : 5000, 10444, 52220000  (nrow, ncol, ncell)
resolution  : 0.01, 0.01  (x, y)
extent      : -26.10815, 78.33185, -10, 40  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Users/hannesmueller/Documents/NDVI Data/MOD13Q1_2000-02-18.250m_16_days_NDVI.tif
names       : MOD13Q1_2000.02.18.250m_16_days_NDVI
values      : -32768, 32767  (min, max)

> r.dwork[[1]]
class       : RasterLayer
band        : 1  (of  1380  bands)
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : /Users/hannesmueller/Dropbox/Locus Project Augustin Lavinia Bruno/gis data/spei/spei01.nc
names       : X1901.01.16
z-value     : 1901-01-16
zvar        : spei

in the band element of the second layer, whose meaning I do not understand. Finally, I exported them with

writeRaster(r.works, "rworks", format = "GTiff")
writeRaster(r.dwork, "rdwork", format = "GTiff")

and the former opens perfectly in a GIS software, whereas the latter does not. That makes me very suspicious about the second raster, which was obtained with the .nc file. Am I importing it into R wrongly and/or missing any step when using zonal statistics/exporting it as geotiffs?


texts: ArcGIS Server Manager (10.2.2) - permissions applied to folder not "sticking"

texts:
I was able to do this before updating to 10.2.2, because I have several other folders with specific permissions applied via local role store.  User store is Active Directory.  Data store is flat files via UNC to shared drive.

Each time I edit the permissions of a newly-created folder, I choose private, un-check allow access to all users, & add the role to allowed roles.  Saving says that, "permissions applied successfully," but the lock icon remains open, & checking the permissions shows that the folder is public.

I have tried deleting & re-creating the folder in the server manager.


texts:
I'm trying to make tool of change detection between two polygons layers by fields of initial state and final state. Everything is working except attribute selection, I would like user could choose whether areals with no change include or exclude from calculations. I create new layer by intersection, then I want to compare values in field of initial state with values in field of final state and exclude areals without change.

# select areals without no change
if change == "Exclude":
    where_clause = '"' + field_initial + '" <> ' + "'" + field_final + "'"
    arcpy.Select_analysis(output_fc, output_select2, where_clause)

output_fc is layer created by intersection.

It creates new layer, but with all objects, it makes no selection at all, I would like to get same result as when I use attribute selection in ArcMAP. I tried to do it in python window, also in tool, there is no error in expression but no selection.


texts: gdal_grid has some problem with reading vrt file

texts:
I am using gdal_grid module(from gdal 1.10) and it has some problem with reading the vrt file.
I checked this link but it does not work. gdal_grid problem reading vrt file

The command is

  ./gdal_grid -l grid grid.vrt newgrid.tif

The Error Message is

  Unable to open input datasource "grid.vrt".

The vrt file is like following

<ORGVRTDataSource>
    <OGRVRTLayer name="grid">
        <SrcDataSource>grid.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="x" y="y" z="z"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

The csv file is like following

x,y,z
565650,5121960,1048
565680,5121960,1043
565710,5121960,1039
565740,5121960,1038
565770,5121960,1036
565800,5121960,1024

Thanks.


texts: Free GIS software for orthorectification?

texts: Could someone recommend a free and/or open source GIS application that can orthorectify map images?


texts: infotemplate getting cutoff on top or bottom of screen?

texts:
in my app, i display few points across the map. these points have infotemplate with details in them and sometimes are quite large (around 350px height) if they have alot of information associated with it. so sometimes on smaller resolutions, either header or footer gets cutoff. is there any way to avoid this?

i am not setting up my anchors, just leaving it up to the arcgis javascript api.

map.infoWindow.show(evt.screenPoint, map.getInfoWindowAnchor(evt.screenPoint));


texts: Strange interpolation with natual-neighbor method

texts:
I have base shapefile and try to use different interpolation methods to get image like this:

This layer is constructed using point to raster conversion.

Generally, interpolation methods as idw and so on works normal, but exactly natural neighbor does very strange thing:

Result is veeeeery long image. What can it be?

The 180th meridian goes through my map, may it be the reason?


texts: Geocoding - looking to identify county from list of zipcodes

texts: Using QGIS, how do I find county of an associated zipcode? I have zipcode for each point, however, I do not have a street address. I am hoping to aggregate my data to the county level using the zipcodes.


texts: Preparing multi class from multi value of raster data?

texts:
How can I prepare a specific band for each value in the below photo?

This raster file is a result of classification which I prepared in ENVI and for some edits I imported it in ArcGIS. after modification I imported the TIFF file  in ENVI for providing confusion matrix, but it was imported as 1 band. I need multi bands of this file to provide a confusion matrix in ENVI.

I am beginner in using ArcGIS and it could be a great help to guide me and give me a solution.


texts: What does Overlay files (*.ovl files) contain about GIS data and how to read them?

texts:
Recently, I started dealing with Overlay files (*.ovl files) and my focus is to read the data from OVL files and display them.

Could you please anyone guide me what are these files and what is the best way to read the content both via product and via programmings (pointers to any knowledge base is very appreciated).

Please let me know if further information is needed.


texts: QGIS: Saved project crashed and I lost three days of work

texts:
I have drawn over 5000 polygons on a new layer, which following a program crash today seem to all be gone. This is despite me saving regularly. I thought that because QGIS references other layer files, my saves would update those layer files with the information I had drawn. Following the crash the respective layer file is 100 bytes and seemingly empty.
If there is any way that I can find the lost information please let me know. I'm running 3.14 on Mac.


texts:
Where can I acquire global wind (speed and direction) data?

Preferably in some kind of GIS format and free of charge.


texts: Save GNSS accuracy for each vertex of polygon/line feature in QField

texts:
In QField it is possible to save positioning information in point features through expression variables like @position_horizontal_accuracy or @position_number_of_used_satellites. The explanation is given in the official QField documentation (https://qfield.org/docs/en/prepare/attributes-form.html).
For a comprehensive data accuracy assessment I want to save the abovementioned variables and maybe additional variables (e.g. @position_vertical_accuracy) for each vertex of line and polygon features. In 2017 Matthias Kuhn mentioned the possibility to theoretically do so by saving the information to the m-value of vertices but that it was not yet implemented (https://github.com/opengisch/QField/issues/179).
Is it possible to somehow store multiple GNSS position information for each vertex of a line/polygon feature? How?
If not, is it at least possible to only store horizontal accuracy as a single value in the m-value of vertices?
If not directly possible, are there any workarounds?


texts: Printable Esri symbol sets, a.k.a. style references

texts:
Does anyone know of good printable documents of the as-shipped Esri symbol sets, also called Style References?

The new with Arcgis 10 ability to search symbol sets is a great boon, but one can still only find things for which you know the keywords for, and only few so many at a time on screen. I find hard copy of things like this much easier and faster to look through and locate the right one for the map at hand.

For example this morning I was looking for a good glacial moraine symbol and searching the style sets didn't locate anything relevant. So I decided to bite the bullet and build my own. In the process I went browsing through the ...\Desktop10\Styles directory and discovered several pictures I've never noticed represented in any styles that would be perfect for the project I'm working on. Since they're here they very likely already exist in a symbol set and I'd be wasting my time creating them anew. Thus the desire for a printed set to help me find them.

In previous versions there is a developer sample called Style Dump, but it doesn't work in Arcgis 10 (and the results were pretty crude and wasted paper anyway).


texts: Avoiding/getting rid of topological errors after dissolve?

texts:
After dissolving a layer in QGIS 2.8, I get a lot of topological errors in the new (dissolved) layer. Any suggestions as to how to get rid of these errors without manually editing each one?

update:

trying to use @Josephs suggestion I get a strange error message:


texts: Set center coordinate for Leaflet-tilelayer-mask layer

texts:
I am working with this Leaflet plugin: Leaflet-tilelayer-mask  and I want to specify the mask center coordinate as opposed to it picking the base map's center coordinate. Tips on how to achieve this?
let map = L.map("map", {
    zoom: 16,
    center: [35.9734, 117.7926]
});

// background layer -- color filter plugin
let myTileLayer = L.tileLayer.colorFilter('https://{s}.tile.thunderforest.com/pioneer/{z}/{x}/{y}.png?apikey=xxxxxxxxxxxxxxxx', {
    attribution: ''
}).addTo(map);

// mask layer -- tileLayer mask plugin
let mask = L.tileLayer.mask('https://{s}.tile.thunderforest.com/pioneer/{z}/{x}/{y}.png?apikey=xxxxxxxxxxxxxxxx', {
    maskSize: 512,
    maxZoom: 18,
    maxNativeZoom: 17
}).addTo(map);

map.on("load", function(e) {
    // mask.setCenter(e.containerPoint);
    mask.setCenter({ lat: 35.9665, lng: 117.8089 }); // this doesn't work

});


texts: Improving QGIS raster calculation performance?

texts:
I'm performing raster analysis on a QGIS 2.14 project CRS is EPSG: 32632.
I used the official 10m DEM supplied by local authorities of Sardinia extracting the area of interest of 446 raster cells.
I've made a slope map in 5 classes, and this task was accomplished without problems.

Then I used Raster Calculator to reclass the slope map assigning weight to each class as follows:

0-10%  -> 2; 11-20% -> 1; 21-35% -> 0; 36-50% -> -1; >50%   -> -2

So I inserted the following expression:

( "Slope_clip_UTM@1" < 11 )  * 2 +  (  ( "Slope_clip_UTM@1" >= 11 )  AND  ( "Slope_clip_UTM@1" < 20 )  )  * 1 +  (  ( "Slope_clip_UTM@1" >= 20 )  AND  ( "Slope_clip_UTM@1" < 36 )  )  * 0 +  (  ( "Slope_clip_UTM@1" >= 36 )  AND  ( "Slope_clip_UTM@1" < 50 )  )  * -1 +  ( "Slope_clip_UTM@1" >= 50 )*-2*

First I've performed raster calculation on a slope raster the CRS of which was EPSG 3003 by default and I obtained an empty raster (no cells and no values), but when I carghed Min/Max values they were both equal 2.

Latter attempt was made right after saving the slope raster in EPSG 32632, the same of the Project: this time as a result I obtained an uniform raster in wich the values of each cell is equal to 0; but when I charged Min/Max this values are both equal to 1.79769e+308.

How do I fix the problem?


texts: Get a large area landsat 8 metadata from earth-engine

texts:
Obviously the code below will returned null because there's no picture in landsat 8 will contain such a large area.

#python
geometry = ee.Geometry.Polygon(
    [[73.0, 25.0],
      [73.0, 49.0],
      [125.0, 25.0],
      [125.0,49.0]])
landsat = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterDate(date_start, data_end).filterBounds(geometry)

But I wonder is there any way that can get a Collection contains pictures which can cover the whole area. I'd prefer these pictures are all metadata before mosaic. Mosaic() make no sense because I want to do some pre-process on single picture.


texts: Leaflet map click event: how to update bar plot in external div?

texts:
I made this code to trigger a graph every time you click on a GeoJSON feature, so that for each polygon in the map, the graph shows data for the corresponding name of the region ("Comune"). The data is contained in a variable called myJSON.

Example:

var myJSON = [{

"": 0,
"Comune": "BONDENO",
"PUNTEGGIOSCUOLA1516": 4.25,
"Value 1": 63,
"Value 2": 8,
"Value 3": 17,
"DANNO": 6,
"Somma valori": 88,
"numero_di_scuole": 4,
"CERTIFICATOAGIBILITAABITABILITA": "",
"DOCUMENTOVALUTAZIONERISCHIO": "",
"PIANOEMERGENZA": "",
"PUNTEGGIOSCUOLA1617": 5.5,
"VINCOLIPAESAGGIO": 1,
"EDIFICIOVETUSTO": 0,
"PROGETTAZIONEANTISISMICA": 1,
"VINCOLIIDROGEOLOGICI": 0,
"PUNTEGGIO_SCUOLE_PARITARIE_1516": "",
"PUNTEGGIO_SCUOLE_PARITARIE_1617": ""]}

When I click on "BONDENO" on the map, the bar plot should show "Value1" and "Value2". However, I can't get it to work.

Here is my code:

function geo_json_dc64d985f99647e7b35676bc82a64cc7_onEachFeature(
  feature,
  layer
) {
  layer.on({
    mouseout: function(e) {
      geo_json_dc64d985f99647e7b35676bc82a64cc7.resetStyle(e.target);
    },
    mouseover: function(e) {
      e.target.setStyle(
        geo_json_dc64d985f99647e7b35676bc82a64cc7_highlighter(e.target.feature)
      );
    },
    click: function(e) {
      map_ae26379126cd4ce79aea9d0f395ec09f.fitBounds(e.target.getBounds());

      myJSON.forEach(function(arrayItem) {
        if (
          arrayItem["Comune"] ==
          e.target.feature.properties.NOME_COM.toUpperCase()
        ) {
          $("#welcome").remove();
          $("#information").html("");
          $("#information").html(
            '

Zona ' +
              arrayItem.Comune +
              '

 punteggioscuole: ' +
              arrayItem.PUNTEGGIOSCUOLA1516
          );
          $("#chartContainer").append(
            '<canvas id="redditiChart" width="200" height="200"><canvas>'
          );
          var grafobase = document
            .getElementById("redditiChart")
            .getContext("2d");

          new Chart(grafobase, {
            type: "bar",
            data: {
              datasets: [
                {
                  label: ?,
                  data: ?,
                  backgroundColor: "rgb(255, 99, 132)",
                  borderWidth: 1
                }
              ]
            }
          });
        }
      });
    }
  });
}

How can I do that? You can see how my code currently works here.


texts: qgis2web leaflet: change position of layers control box

texts:
I exported a leaflet-webmap with qgis2web. On the top-left, there are the controls (like zooming, serching etc.). On the top-right, there is the layers control box (checking/unchecking the layers of the map). This is the one I want to change.

I want to move the layers control box to the left side, under the other control-box. When checking the source-code in my browser, the index.html shows a line <div class="leaflet-top leaflet-right"> - however, if I look to my saved index.html, there is no such line - it seems to be generated from somewhere else. But from where? I checked the files in css and js folders - found nothing. How can I change the position of this box?


texts: QGIS: sip module implements API v11.0 to v11.3 but the qgis._core module requires API v10.0

texts:
I have read all the questions on the SIP module and QGIS, but they do not answer my question.  Most of these are dealing with screwed up paths; I do not think my environment variables are incorrect.

I've got a Windows 7 x64 system running Boundless Desktop 1.0 (QGIS 2.14 Essen).  I set its python install to be the primary system python, and have been pip-ing a bunch of libraries I need.  Here are my user environment variables:

GDAL_DATA = C:\Boundless\osgeo4w\share\gdal
GDAL_DRIVER_PATH = C:\Boundless\osgeo4w\bin
PATH = [path variables],C:\Boundless\osgeo4w\bin;C:\Boundless\osgeo4w\apps\qgis\bin;C:\Boundless\osgeo4w\apps\Python27\Scripts
PYTHONHOME = C:\Boundless\osgeo4w\apps\Python27
PYTHONPATH = C:\Boundless\osgeo4w\apps\Python27;C:\Boundless\osgeo4w\apps\Python27\Scripts;C:\Boundless\osgeo4w\apps\Python27\DLLs;C:\Boundless\osgeo4w\apps\Python27\Lib

When I start QGIS, I am faced with this:

RuntimeError: the sip module implements API v11.0 to v11.3 but the qgis._core module requires API v10.0

When I check sip within a python shell, I get the following:

>>>import sip
>>>print(sip, sip.SIP_VERSION_STR)
(module 'sip' from 'C:\Boundless\osgeo4w\apps\Python27\lib\site-packages\sip.pyd'>, '4.18.1')

When I check sip from the command line, I get the following:

C:\>sip -V
4.14.7

So I've got a mismatch between two versions of sip.  The QGIS error indicates that it is looking for an older version.

I have, since installing and configuring QGIS, updated GDAL (to 2.1.3) and PyQt4 (to 4.8.6).  My feeling is that these actions have broken QGIS, and my question is this:

Is there a particular version of GDAL/PyQt4 I can just downgrade to, or would the better option be to reinstall QGIS, thereby (hopefully) overwriting and replacing the offending libraries?


texts: Assigning row number for sorted table by column in QGIS

texts:
Field calculator has a great function under "Records" called "$rownum" which returns (assigns) the number of the current row.
Is there a way to make it work together with sorting?
Let's say I want to sort table by column "length" and get the order (sequence) of lines by length from shortest to longest.
I know I can use some office spreadsheet software easily for this to sort by column and populate new column with order (sequence). But I would prefer to do the whole job in QGIS.
As workaround there is also Changing order of features in shapefile - using mmqgis plugin and function modify/sort and after that use "$rownum" on the new saved file.


texts:
I used gdal_translate to converted a Float32 DEM to a UInt16 .tif and then ran gdalinfo -hist. In the histogram it says that it has 256 buckets. I also ran -hist on a 8bit .tif and that also has 256 buckets. Does that mean that the values doesn't span over the whole spectrum in my 16bit .tif and if so what would I do to make more of the range used?
The min and max in the histogram is 0 and 65535.

I want to get as much values as possible in my 16bit as I plan to use it as a heightmap/displacement map later and want as much of the original detail in there as possible.

This is what I ran with gdal_translate

gdal_translate -of GTiff -ot UInt16 -scale -5660.574 -5618.043 0 65535 input.tif output.tif


texts: How to geocode US addresses using gisgraphy (400k)

texts: Would someone please direct me in geocoding around 400k addresses using gisgraphy or openstreetmap?


texts: How can I increase ArcMap font size, I can barely read it!

texts: I'm sure there is a simple answer out there, but how do I increase the font size in the ArcMap work windows?  I am using Windows 8 on a Surface tablet.


texts:
I'm curious to know what the correct projected CRS for Uganda is and if there is a national grid?
Is it ok to use UTM as it is conformal?


texts: Where are the GeoServer HTTP logs on Jetty server?

texts:
I am trying to set up a load balanced cluster of GeoServers using Google Cloud, this requires setting a "health check". I'm trying to health check by doing an http get on /geoserver/web/ to load the main admin page for no other reason to see that the server is alive. Apparently this times out and fails. In trying to debug this, I can manually load [server:port]/geoserver/web but on the server side I do not see any HTTP access logs to see my manual access or the health check attempt access. I launch my server with the command

java [LIST OF OPTIONS] start.jar > log.txt

I can see startup information in log.txt and when I do WMS fetches those show up there too. But I don't see regular http accesses on the admin page. Where can I find/log those?


texts: Using ogr2ogr to convert all shapefiles in directory?

texts:
I have a directory with several shapefiles.

How can convert all these shapefiles to MapInfo with ogr2ogr?

I know how I can convert one file.
And I can make a batch script writing a line for each file.
But isn't there an easier way to convert all the files in a directory (and subdirectory).


texts: Transforming coordinates from WGS84 to Gauss-Krüger 5

texts:
GPS coordinates in the format of WGS84 are quite common. I need to transform them into Gauss-Krüger coordinates but I have no idea how to do that. The only things I know are that it is not easy and I have to use the Helmert-transformation.

I usually work with r-project and would therefore prefer to work with a r-package to solve my problem. Is there a package which could handle that?

Is it actually necessary to transform the "GPS coordinates" in order to create a hillshade map with different "landmarks" according to the coordinates? Is there any open source software to create a map and then mark the positions based on the coordinates?


texts: Changing format of column from character to integer using MapInfo?

texts:
I have two tables which I am trying to do distance calculations on in the SQL select window of MapInfo. I have completed these calculations before so I am happy that the formula I am using is correct, however; MapInfo shows a mismatch error which I understand is to do with the colums I am trying to match not being the same format i.e. one student ID is in string (Character) and one student ID is in value (integer). I have checked this in the table structure and this is the case. I have tried to change the format in the CSV and overwrite the table but the issue remains. I have also tried adding: int(table.stud_id) to the WHERE clause and this clears the "data mis-match" error but I then get an "invalid join in WHERE clause error".

Is there away to change the format of the column in the table from character to interger in Mapinfo?


texts: Determining Thalweg line of river without bathymetry?

texts:
I have the boundaries and centerline of a river.

The Thalweg line defines the deepest elevation of a watercourse (see Image taken from Google). I do not have accurate enough  bathymetry or Digital Elevation Models to determine the proper deepest elevation of a river channel.

I would like to know if it is possible determine the Thalweg line of the river without the bathmetry?

It looks like the Thalweg line follows a predictable pattern where it hugs the outside curve of a river when it meanders and is in the center of the river with it is straight.

I am using ArcGIS 10.1 with all licenses available, and I am familiar with python.


texts: Calculate percentages with field calculator in arcgis?

texts:
I am trying to calculate the percentage of each the area of each polygon in the selected area.

I tried using ( math.fsum(!shape.area!) / !PolygonArea! )*100 and other codes but I cannot make it work.

Suggestions on how to solve this?


texts: Understanding why ArcPy Cost Path Analysis tool faster than ArcObjects?

texts:
Although I use python to create geoprocessing scripts/services, I was under the impression that using ArcObjects to do the equivalent operation(s) will have better performance.

I've posted ArcGIS Server GP Service - RasterIO.dll crashing ArcSOC.exe and ArcGIS Geoprocessing Script runs fine in Desktop but crashes as Geoprocessing Service? over the past couple of days about getting geoprocessing scripts that use Spatial Analyst tools to work as geoprocessing services.  My deadline is fast approaching, so I have decided to go the SOE route to achieve the desired functionality.

Getting a cost path analysis in ArcObjects was relatively straight-forward using the .NET ESRI.ArcGIS.SpatialAnalyst.RasterDistanceOpClass, specifically the CostDistanceFull() and CostPath() Methods.

Some code snippets of how I am doing things:

  Python

# Get Cost Path Origin and Destination Points
inputPointsShp = 'D:/RasterStuff/test_points.shp'
arcpy.MakeFeatureLayer_management(inputPointsShp,"origin",' "TYPE" = \'ORIGIN\' ')
arcpy.MakeFeatureLayer_management(inputPointsShp,"destination",' "TYPE" = \'DESTINATION\' ')

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute CostDistance
outCostDistance = CostDistance("origin",SOURCE_RASTER,"#","backlink")

# Execute CostPath
outCostPath = CostPath("destination", outCostDistance,"backlink")

# Convert Result to Polyline
arcpy.RasterToPolyline_conversion(outCostPath, "leastCostPath")
featSet = arcpy.FeatureSet("leastCostPath")

  C#

IDistanceOp distanceOp = new RasterDistanceOpClass();
IRasterBandCollection costDistanceRaster = (IRasterBandCollection)distanceOp.CostDistanceFull((IGeoDataset)sourceFc, (IGeoDataset)raster, true, true, false);
IRasterBand distanceRaster = costDistanceRaster.Item(0);
IRasterBand backLinkRaster = costDistanceRaster.Item(1);

IGeoDataset costPath = distanceOp.CostPath((IGeoDataset)destFc, (IGeoDataset)distanceRaster, (IGeoDataset)backLinkRaster, ESRI.ArcGIS.SpatialAnalyst.esriGeoAnalysisPathEnum.esriGeoAnalysisPathForEachCell);

A cost path analysis in ArcPy (using sa.CostDistance and sa.CostPath) takes approx 15-20 sec.  Using the exact same inputs, the ArcObjects based routine takes 55-60 sec.  Even using the .NET Geoprocessor is significantly slower than arcpy.

I guess my questions here are:

Are the ArcPy and ArcObjects implementations pointing at the same code base (through their Python and .NET wrappers)?
Any tips to optimize ArcObject based Cost Path analysis?


texts: How to indicate flow of a watershed feature?

texts:
While tracing and digitizing rivers, streams, rivulets, etc., is it possible to indicate the direction/flow in the shapefile on QGIS?

I'd like to avoid using map labels. Ideally, when done tracing, I would like to change the symbology to a line pattern of tiny arrows showing the flow direction.

Some similar questions I've found here deal with elevation and contours obtained from DEM materials. Although I am also working with DEM files, I would like to be able to work this without using the ASTER material and simply input it myself since I know the flow direction.


texts:
There are 4 required inputs for this addon (r.damflood) in GRASS.

elev=  elevation raster map (raster file with elevation data)
lake=  water depth raster map (raster file of only the lake area)
dambreak= Name of dam breach width raster map
manning= Name of Manning's roughness coefficient raster map

For the last two (dambreak and manning), I have had some trouble when trying to get them.
First, for the dambreak file, the idea is to do 'breach = DTM - your map'. But how can I do this with only one raster file?
Second of all, for the manning map, I am supposed to use the r.reclass tool. But what could I possibly classify in an elevation map?


texts: Test Attribute with 1 or more possible values in FME?

texts:
This question continues story on link What's the syntax for an AttributeFilter in FME?.

Now, how can I test attribute value if there are more than one value and proceed it to the Clipper (uses as a multiple clippers)? Use some looping method?

In xml file is for example

So there can be maybe only one or can be many of them.

I tried and Tester output is only one value (last of possible values, in my example above Tester output is only areaID=1253).

Tnx,


texts:
There is a PostGIS database that has a road network imported from Open Street Map, stored in osm_roads table.

For the testing purposes I need to change some road geometry there. I tried to edit it with uDig Version: 1.4.0b. When I select the layer with osm_roads, I don't see any 'edit' buttons; the only buttons I see are 'create arc', but after drawing the arc it disappears, so I can't edit the road network anyway.

Please advise what I should configure to make this changes possible. What is the most likely reason why uDig doesn't allow to edit the road layer?


texts: Inserting script into ArcGIS ModelBuilder?

texts:
I am looking for some tips about inserting a script in the ModelBuilder...

I am actually building a model in the ArcGIS ModelBuilder. At one point I would like to insert a script that crosses information from two layers at one point of the model with a output layer, that will be used further along.

In ArcMap, I added the script with a right-click in My Toolbox > Add Script...

After loading my script, I applied the following parameters:

Two inputs, referring my layers that I use in the script (raster layer type)
One output, for the result (raster layer type)

I then go to my model and drag my script into it. I connect my two layers as inputs with the arrow and specify the output file name. All seems ok as the boxes appear with the correct colors.

I have to specify that I am working with a variable folder %folder% so I added that in my script.

When I validate my model it goes all blank, because it cannot find the input layers, as they are created in the model.

Can anybody help me please?

Here is my script:

# ---------------------------------------------------------------------------
# calculator.py
# Created on: 2012-11-15 16:59:11.00000
#   (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Local variables:
raster = "%Dossier%\Erosion.gdb\raster"

# Process: Calculatrice raster
arcpy.gp.RasterCalculator_sa("Con((\"%Surf_Drainee%\" > 12) & (\"%Surf_Drainee%\" <= 80) & (\"%Pente_Rclass%\" <= 24), \"%Pente_Rclass%\" + 1, Con((\"%Surf_Drainee%\" > 80) & (\"%Pente_Rclass%\" < 24), \"%Pente_Rclass%\" + 2, Con((\"%Surf_Drainee%\" > 80) & (\"%Pente_Rclass%\" == 24), \"%Pente_Rclass%\" + 1, \"%Pente_Rclass%\")))", raster)


texts: Is it possible to add an overlay from Google Earth to QGIS?

texts:
For a hiking project in Kyrgyzstan we're planning our trails in GoogleEarth. We georeferenced .gif-files into an overlay-file (.kmz-file).
How can I import this file into QGIS?
If I us "Add Vector", I only get a shapefile, where the maps are shown as monochrome squares.


texts: Calculate focal statistics from a 10m raster to populate a 1km raster, using the 1km raster to define the window?

texts:
I have a binary raster (1-0), cell size = 10m.
I would like to generate a 1km raster, where cell values equal the mean value of all cells within a 3 km radius. I know I can use focal statistics to generate a 1km raster from a 1km raster using a window.
But is there any way to calculate the mean of all the 10m cells that fall within a window around each 1km cell (i.e. window being a circular with radius = 3km)?


texts: Adding intersecting information from point layer to polygon layer pop-up window in ArcGIS Online using Arcade

texts:
I have two layers. One point layer and one polygon layer.
The point layer has several points that fall within each polygon. The points represent different types of sports, and the polygons represent the buildings that house these sports.
I want to be able to click on a building, and have a pop-up which lists all the available sports within that building. I am using Arcade expressions in the pop-up, but am stuck on only being able to list one of the sports on the pop-up using this tip:
https://community.esri.com/community/gis/web-gis/arcgisonline/blog/2018/12/10/overlapping-features-in-pop-ups-quick-introduction-to-using-featuresets-with-arcade
Is there a way to list all of the sports that intersect?


texts:
I have created a blank print layout, now I want to add items to it like a map, legend, etc.

Here is my code so far to create a new print layout, this works just fine:

#This creates a new print layout
project = QgsProject.instance()
manager = project.layoutManager()
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName('My Layout')
manager.addLayout(layout)

Now I want to add stuff to the layout, starting with a map:

map_item = QgsLayoutItemMap.create(layout)

But I get the following error:

When I look at the documentation for QgsLayoutItemMap.create I see this, except I do not understand it:

What does this mean "the caller takes responsibility for deleting the returned object" and how do I add a map (or other) objects to my Print Layout?


texts: FME - Use Python Caller to Set Search Envelope in Another Workspace

texts:
I have an FME workspace which has a number of different workspaceRunner transformers.  The first workspaceRunner launches a workspace which buffers a pipeline and outputs the bounding coordinates to a CSV file.

I then have a python caller which reads the CSV to get the bounding coordinates and save them as global variables.

The code of the Python Script looks like this:

import fmeobjects
import csv

# Template Function interface:
def processFeature(feature):
    global XMIN
    global XMAX
    global YMIN
    global YMAX

    outputDirectory = FME_MacroValues['Output_Directory']
    csvPath = outputDirectory + '\\Buffer Bounds.csv'
    coordlist = []

    with open(csvPath, 'rb') as csvfile:
        reader = csv.reader(csvfile, delimiter = ',')
        for row in reader:
            coordlist.append(row)

    XMIN = float(coordlist[0][0])
    XMAX = float(coordlist[0][2])
    YMIN = float(coordlist[0][3])
    YMAX = float(coordlist[0][4])

    print XMIN, XMAX, YMIN, YMAX
    pass

# Template Class Interface:
class FeatureProcessor(object):
    def __init__(self):
        pass
    def input(self,feature):
        self.pyoutput(feature)
    def close(self):
        pass

After the PythonCaller, I have four VariableRetriever transformers which get the values of XMIN, XMAX, YMIN, YMAX and save them as attributes.

That is then connected to a workspaceRunner transformer.  The parameters for the bounding box are passed to the next workspace and are used to define the Search Envelope Coordinates for the File Geodatabase.

My problem is that these variables are somehow not getting passed to the attributes.  Am I using the variableRetriever transformers correctly?  I have read on some FME forums that you can use them this way.  I've declared them as global variables.


texts: Seeking automated approach for categorising roads according to their designs using QGIS?

texts:
I have a shapefile comprised of polylines. The road segments are represented by lines that are between nodes. I want to categorise each road segment accordingly as they appear either as a straight road, curvilinear road or a road with single (or multiple) blind turning(s).

I want to refrain from using a visual approach when categorising these segments because its labour intensive, and such approach is subjective which may introduce some form of misclassification bias.

I was wondering if there is an automated approach for performing this task in QGIS?


texts: Seeking open source alternative to RivEX (add-in for ArcGIS designed to process vector river networks)

texts:
I wish to process a river network - for the assessment of hydrological potential. RivEX is one addin (not free) which could be used with ArcGIS.

Is there any alternative to this?


texts:
I am developing a little web app to transform coordinates between CRSs. To do that I am using proj4js, but I notice that proj4js can't work with geocentric CRSs.

Is there another tool similar to proj4 to work with geocentric CRSs, or I have to code the formulas between geocentric-geographic coordinates. Once I have geographic coordinates I use proj4js to obtain projected coordinates.


texts: Geoserver: Georectification via SLD fails

texts:
Situation

I have a raster image and the coordinates of its four corner pixels. My goal is to warp the image with respect to this coordinates ("georectify") on the fly when requested via GeoServer to display in OpenLayers.

Inspired by this discussion list, my idea was to do the following in GeoServer:

Create a WorldImage layer (raster + text file with geodata)
Create a Styled Layer Descriptor (SLD) calling the Web Processing Service (WPS) GeorectifyCoverage which "georectifies a raster via Ground Control Points using gdal_warp".
Apply the SLD to the layer.

Problem

When I request the layer for example via the OpenLayers preview, it throws an error reporting:

Could not convert the value EPSG:32632 into the expected type interface
org.opengis.referencing.crs.CoordinateReferenceSystem for parameter targetCRS

Oddly enough, when I use the "WPS request builder" (in the GeoServer Web Application) to call GeorectifyCoverage with the same parameters, it gives the right output.
So the error must originate somewhere from my SLD document below.

Can anybody give me a clue what might be wrong with it?

Additional Info

SLD document calling WPS rendering transformation GeorectifyCoverage:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
    <NamedLayer>
        <Name>GeorectifyCoverage</Name>
        <UserStyle>
            <Title>Georectify Coverage</Title>
            <Abstract>Warps image to specified GCP</Abstract>
            <FeatureTypeStyle>
                <Transformation>
                    <ogc:Function name="gs:GeorectifyCoverage">
                        <ogc:Function name="parameter">
                            <ogc:Literal>data</ogc:Literal>
                        </ogc:Function>
                        <ogc:Function name="parameter">
                            <ogc:Literal>gcp</ogc:Literal>
                            <ogc:Literal>[[[0, 0], [308639.500, 5154180.500]], [[0, 450], [308639.500, 5150929.500]], [[600, 0], [315730.500, 5154180.500]], [[600, 450], [315730.500, 5150929.500]]]</ogc:Literal>
                        </ogc:Function>
                        <ogc:Function name="parameter">
                            <ogc:Literal>targetCRS</ogc:Literal>
                            <ogc:Literal>EPSG:32632</ogc:Literal>
                        </ogc:Function>
                    </ogc:Function>
                </Transformation>
                <Rule>
                    <RasterSymbolizer>
                        <Opacity>1.0</Opacity>
                    </RasterSymbolizer>
                </Rule>
            </FeatureTypeStyle>
        </UserStyle>
    </NamedLayer>
</StyledLayerDescriptor>

where GeorectifyCoverage requires the following input:

data: input raster (my layer name)
gcp: List of Ground control points (my corner coordinates)
targetCRS: system to use for the output raster (my CRS corresponding to the GCP)

Notes:

When I call a different WPS function like Contour using the same document structure and raster image it works perfectly.
When I change EPSG:32632 to EPSG:4326 the same error is reported.
Changing the raster format (jpeg, TIFF) does not helo.
Processing the same data offline with the command line tools gdal_translate (assign GCPs) and gdalwarp results in the desired output.

My System:

Ubuntu 12.04
OpenGeo Suite 3.0.2
GeoServer 2.2

Alternative approaches using OpenGeo would also be welcome.


texts: Geosever : make line string feature from MongoDB collection

texts:
I installed the gt-mongo plugin to a geoserver and made a new MongoDB store and here is the collection documents example:

{
    "_id" : ObjectId("5a2543fcccc8fa54bc0822ac"),
    "loc" : {
        "type" : "Point",
        "coordinates" : [
            22.697815,
            75.740113
        ]
    }
},
{
    "_id" : ObjectId("5a2543fcccc8fa54bc0823ca"),
    "loc" : {
        "type" : "Point",
        "coordinates" : [
            25.697815,
            72.740113
        ]
    }
}

and the schema_store JSON file:

{
    "typeName": "test",
    "userData": {
        "collection": "test"
    },
    "geometryDescriptor": {
        "localName": "location",
        "crs": {
             "properties": {
                 "name": "urn:ogc:def:crs:EPSG:4326"
              },
             "type": "name"
        }
    },
    "attributeDescriptors": [
        {
            "localName": "location",
            "type": {
                "binding": "com.vividsolutions.jts.geom.Point"
            },
            "userData": {
                "encoding": "GeoJSON",
                "mapping": "loc"
            }
        }
    ]
}

I need to serve it as LineString, not as points. any helps?


texts:
I've been trying to run r.sim.water on a 3 meter LiDAR-generated DEM I have, but keep getting G_malloc errors:

r.sim.water --verbose elevin=meterDEM.final@Sediment_Simulation dxin=centDEM.dx@Sediment_Simulation dyin=centDEM.dy@Sediment_Simulation depth=sim.water.depth disch=sim.water.discharge err=sim.water.error nwalk=1000

My error message:

ERROR: G_malloc: unable to allocate 43608 bytes at main.c:488

I'm running GRASS 6.4.2 on Windows.  I've tried using a Ubuntu install on VirtualBox, and it fails out without any error at all.  I tried ratcheting down the number of walkers (hence nwalk=1000), but that doesn't seem to help.  I can't really think of what to try next.


texts:
I have created a IDW heatmap using QGIS 1.8.0 and need to import the raster into PostGIS 2.0.1. QGIS saves the raster layer as a .asc file.

How do you import the raster into PostGIS as a Raster? What format should QGIS save the raster as? Can pgAdmin's PostGIS Shapefile and DBS Loader be used?


texts:
Can the error log created by validate topology be assessed through arcpy? I can validate the following topology -

arcpy.ValidateTopology_management(r"C:\WorkSpace\CWD.gdb\FD\FD_Topology")

Can I access the log created without entering ArcMap? I don't need to be able to fix errors. I just want to know that there are not errors before exporting features. If they exist and I can verify that, I can raise an exception.


texts: pulling an excel sheet with =rand(between) in a Python while loop, and exporting results as .dbf

texts:
To preface my question: I am very new to stack overflow, and relatively new to Python.

I am working on setting up a sensitivity analysis. I am working with 40 parameters that range from 0.1 - 1. My analysis requires simultaneously varying these parameters by +-0.1 roughly ~500 times. These values will then be fed into an ArcGIS tool. So, I need 500 sets of random values, for 40 parameters. These parameter values will then be compared to the output of the tool, to see which parameters the model is the most sensitive to. I've set up an Excel sheet that will randomly calculate these values each time it's opened, but the issue is that they need to be in .dbf format to be read by the ArcGIS tool.

I have set up a while loop (for 10 iterations to start, but will need to be ~500) and tried two different methods, in hopes that I could automate the process of calling the .xls to generate random numbers, and then exporting it to .dbf.

The first, arcpy.CopyRows_management correctly exported to .dbf. The issue was that the output was the exact same for each iteration, and instead of having values of 0.1, 0.2, 0.3 etc, it contained values of 0.22, 0.37, 0.68 etc. It wasn't to the tenths, even though that was specified in the formulas in the .xls.

I also tried arcpy.TabletoTable_conversion but that was throwing the ERROR 999999: Error executing function.

I am open to all kinds of suggestions. Perhaps there is an easier way to randomly sample and export results to .dbf in Python. This does not need to be done using arcpy, but that is all I've really worked with. I really appreciate any help that is provided! Thanks for your time.

i = 0
while i < 10:
    # Set run specific variables
    lulc = "D:\\SARuns\\lulc_nosoils_rand.xls\\lulc_nosoils$"
    folder = "D:\\SARuns"

    print "Reading lulc"

    newlulc = "D:\\SARuns\\lulc_nosoils_rand"+str(i)+".dbf"

    print "Reading newlulc"

    # Copy rows is copying it to a dbf, but the values inside
    # are the same for each run. And, none are correct.
    arcpy.CopyRows_management(lulc, newlulc)

    # Table to table should work. But isn't.
    # arcpy.TableToTable_conversion(lulc, folder, newlulc)

    print "Converting table"
    i+= 1


texts:
I have been looking for paleoclimatology data.  I have found some with NOAA ( http://www.ncdc.noaa.gov/paleo/data.html ) but all the data is stored as a KMZ and a series of linking text and html files.  Does anyone know of vendors or sites that have this type of data in a better spatial format?  I am most interested in data around the Great Lakes.

Thanks...


texts:
I want to use the Identify Task with the ArcGIS JavaScript API, similar to this example. I am getting the following errors when I try to add my map service to the code.

  Uncaught TypeError: Cannot read property 'lods' of null

and

  Uncaught TypeError: Cannot read property 'pause' of null

I parsed down the code I am using to what is below.

<!DOCTYPE html>

<meta charset="utf-8">
<meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
<title>Identify Task - 4.3</title>

<link rel="stylesheet" href="https://js.arcgis.com/4.3/esri/css/main.css">

Any ideas on how I can get the map service to display and not throw the errors?


texts: Sum fields from ListFields - RuntimeError: A column was specified that does not exist?

texts:
Does anyone have a suggestion on how to otherwise solve the problem ...
The need to automatically generate all the columns from the table, the calculation of the sum, the results are updated in the new column? Recognizing the conditions: 1) did not know the names of columns 2) the name of the column to update the known
example table:

PS
I am an absolute beginner in python programming

 import arcpy, os, string
 fc = "D:\\path\\test.shp"

 fieldNameList = []
 fields = arcpy.ListFields(fc, "*")

    for field in fields:
       if field.type in ("Double", "Integer", "Single", "SmallInteger"):

          fieldNameList.append(field.name)
          fn = "'" + "', '".join(fieldNameList) + "'"
 print(fn)
 with arcpy.da.UpdateCursor(fc, [fn]) as cursor:

 for row in cursor:

            row[4] = row[0] + row[1] + row[2] + row[3]
            cursor.updateRow(row)

ERROR - RESULT

 >>>
 'tt', 'ff', 'vv', 'rr', 'update'

 Traceback (most recent call last):
   File "E:\path\test.py", line 17, in <module>
   for row in cursor:
 RuntimeError: A column was specified that does not exist.
>>>


texts:
I have a geopandas DataFrame (created from a shapefile) which represents the course of a public transport vessel.

Exemplary DataFrame:

| oid | geometry            |
|-----|---------------------|
|  1  | LineString(1 1, 1 2)|
|  2  | LineString(2 3, 3 4)|
|  3  | LineString(1 2, 2 3)|
|  4  | LineString(3 4, 5 6)|

The DataFrame contains segments of LineStrings in no given order of the object id oid. How can I sort the LineString segments by location to have the joined course segment by segment in geographic order (from example data I want to get the order 1 -> 3 -> 2 -> 4)?

I other words, I want to find the touchpoints row-wise.

Background: I have another DataFrame where I have a measured variable given only by the distance from the starting point of the course (of course along this course) and no GIS information.


texts: update/overwrite feature layer in arcgis.com

texts:
I'm using arcgis.com to grant access to non-GIS staff.

Using the original table data they provided I published a service with multiple feature layers and spent a lot of time configuring their symbology and behavior in the map viewer.

A week later I was informed that they wanted to add a column to the data and that new records would be added to the original data set as time goes on.

Will I have to recreate everything in ArcMap and republish the service and create a new map feature service or is there a way I can overwrite an existing service so that the feature layer can have new features and attribute table updates but keep the symbology in the current map viewer?


texts: How to process large shapefiles faster in ArcGIS for Desktop?

texts:
I have huge line and polygon shapefiles that I want to edit with some geoprocessing tasks (clip, buffer, erase,...). Altough my computer is not that slow, I am annoyed of waiting up to 15 hours for one of these tasks.
Some people say that importing the shapefiles in a File Geodatabase will help. Are there any further hints? Is there any difference whether I perform such ArcGIS for Desktop tasks through the ArcCatalog or ArcMap applications?


texts: Using the ee.Reducer.sum() to count the area of different classes of a layer (hansen/global_forest_change_2017_v1_5)

texts:
I have used the following function to calculate the area of an area of interest (AOI) of a classification using landsat8 images.
This counts the area in square meters of the classification for a given polygon.

var areaChart = ui.Chart.image.byClass({
  image: ee.Image.pixelArea().addBands(lista_imagenes[0]),
  classBand: 'classification',
  region: geometry,
  scale: 30,
  reducer: ee.Reducer.sum()
});

The problem is that the band of classification of landsat8 is just one (classification band). And the layer I am using has three layers (it is call year of loss).

There is any way to count the classes of the layer given a polygon using a similar question?

Here the code:
https://code.earthengine.google.com/495ed315756b7181b4570546be82e2f4


texts:
I'm looking for an efficient way to cluster lines independent of their direction. That means that a line between New York and Los Angeles should be in the same cluster as a line in the other direction between Los Angeles and New York. The start / end point locations should be similar (i.e. San Diego to Long Island should be in the same cluster as LA-NY but probably not San Francisco to Boston) and there are no intermediate points. Input data would be similar to this example:

(By Cassiopeia sweet at Japanese Wikipedia GFDL or CC-BY-SA-3.0, via Wikimedia Commons)

I have previously tried to sort the lines in advance, e.g. to make them all run from west to east, but this does not solve the problem for lines running from north to south and other way around.

Do you know any algorithm dealing with this problem? I've been looking but besides Algorithm to compute average direction of undirected segments I haven't found anything remotely helpful, so I must be using the wrong search terms.


texts: Why can't I translate GeoServer layers to KML?

texts:
I have installed GeoServer 2.16.1 on two independent servers; but both are unable to translate layers to KML (Keyhole Mark-up Language). I receive the following content instead of the KML:

<ServiceExceptionReport version="1.1.1">
    <ServiceException code="InvalidSRS">
        Error occurred decoding the espg code EPSG:2908/kml?layers=nyc:nyc_buildings No code "EPSG:2908/KML?LAYERS=NYC:NYC_BUILDINGS" from authority "European Petroleum Survey Group" found for object of type "IdentifiedObject".
</ServiceException>
</ServiceExceptionReport>

Apparently, the server is trying to interpret "EPSG:2908/KML?LAYERS=NYC:NYC_BUILDINGS" as the EPSG code, instead of just "EPSG:2908".

The URL to generate the KML is:

  localhost:8080/geoserver/nyc/wms?service=WMS&version=1.1.0&request=GetMap&layers=nyc%3Anyc_buildings&bbox=983797.5%2C207443.0%2C991899.0625%2C218850.828125&width=545&height=768&srs=EPSG%3A2908/kml?layers=nyc:nyc_buildings

Any ideas on why this is happening?


texts: Stylized map workflow: Illustrator to Google Maps / OpenLayers

texts:
The aim My current aim is to create a series of (heavily) stylized diagrams and maps. I'm not very familiar with any GIS software, so my weapon of choice would be to use illustrator for creating the maps (the source would be from a GIS system, so the input would be among the lines of a shapefile and is negiotiable).

The problem So far so good, and I do not forsee any problems in this. The main question: how can the workflow from illustrator to any kind of tileserver with the goal to be displayed in google maps / Openlayers or something alike be achieved as painless as possible?

Vector vs Pixels: My idea now is to use a tileserver, because I think rasterizing the images will be easier then to import the styles from illustrator back to a GIS server. However, vectorized data might mean a lot less data to transfer and a better user experience.

Also: Would you recommend this workflow? Am I going about this the wrong way?

Your ideas and suggestions are highly appreciated.


texts: What's the purpose of the reference latitude in Lambert Conformal Conic?

texts:
Going through flight school I got acquainted with the Lambert Conformal Conic projection, as it's the common projection for VFR sectional charts, and I think I have pretty good grasp of it, geometrically.

However, trying to wrap my head around the mathematics behind it I browse around on Wikipedia and on MathWorld which both more or less simply list the formulas, not really diving into the mathematics behind it.

It appears to have 4 parameters:

The two standard parallels,
reference longitude, and reference latitude

What is the purpose of the reference longitude/latitude? Do they affect the actual projection, or will they simply shift the result to let the reference coordinates result in (x,y)=(0,0)?


texts: What is the relationship between the bearing the 'road' distance between two points?

texts:
Given the two points, we know how to calculate the shortest distance between them, and the bearing from one point, in order to reach the other point.

Also, we know the loxodromic distance between them, and the corresponding bearing. In this article, the calculation is explained.

Now I would like to know how to calculate the distance when the bearing is given. This hypothesis means also that you have to give the "straight" road given the bearing, we can suppose that the road is obtained by cutting the sphere with a plan, on the two points. If I am not wrong, you have a unique cut, given the bearing and the two points.


texts:
We are building a spatial application targeted to laypeople which will involve viewing and editing vector geodata (currently just points, but linestrings and polygons in the future) with a raster underlay. We'll deploy both web and mobile (iPhone now, Android and perhaps others later) versions. The server is a custom app we are writing, and there are no ESRI products on the server side.

My question, therefore, is: How do we choose which client-side mapping toolkit to use? I'm aware of the following options, are there others?

Google Maps API or similar (Yahoo, MapQuest, etc.)
OpenLayers (does it work on mobile browsers?)
ESRI's web/mobile APIs

We are particularly interested in the pros/cons of the ESRI APIs under this scenario. Our company does have some kind of relationship with ESRI, though I know nothing about the details (and our group currently has no involvement with ESRI or ESRI products).


texts: How to get a layer object from a geometry

texts:
I have a polygon as a gdal geometry and I'd like to use it with  gdal.RasterizeLayer(). As its name states, this method needs a layer object, not a geometry. So I need to convert/wrap my geometry into a layer object.

The only way I got this to work is by exporting the geometry as GeoJSON and opening that with the GeoJSON driver as datasource:

GeoJSONdriver = ogr.GetDriverByName('GeoJSON')
ds = GeoJSONdriver.Open(my_polygon_geometry.ExportToJson())

layer = ds.GetLayer(0)

This feels very clumsy and since I also need to look out for performance I was wondering if there is a more direct approach to make that conversion?


texts: Toggle between Layers in QGIS Layer Panel

texts:
Working on city mapping...
We've several layers that are constantly being edited and saved.

Would like to be able to save a few clicks by creating a custom macro/hotkey to "save current layer", and scroll up or down to the next layer in the Layers Panel.

Thinking "ctrl+mouse wheel up or down".

So far have been unable to find anything of this sort.


texts: Can we create relationship classes based on a geometric network?

texts:
I have a point Feature Class and a polyline Feature class which I connected in a geometric network. I also need to connect them logically. In my theoretical model, an edge is defined by the two points at its end. So I need a Many-to Many relationship. An edge can have two points and a point can have two edges. Connecting them in a relationship requires the maintenance of unique ID's etc.

So I was wondering is there a possibility that based on their geographical relationship of points and line to create as well a logical relationship in the database, and save me from the whole process of giving them ID's, connecting them, etc?

It's probably a rather general and abstract question, but I guess some of you might understand what I mean.

Thank you!

I use ArcGIS 10.0.


texts: Style line-width proportionally to maintain relative size in Mapbox GL

texts:
I want a data driven line-width to keep its relative size in all sizes
I tried using an exponential function for that but it did not work the way I expected.
For the sake of the example let's say I have 3 lines with value: 1, 10, 100

In zoom 0 their width would be: 1px, 10px, 100px
In zoom 1: 2px, 20px, 200px
In zoom 2: 4px, 40px, 400px
In zoom 4: 16px, 160px, 1600px
In zoom 10: 1024px, 10240px, 102400px
etc…

I marked the rate of change to 1.99 which gave the closest results, but it's still not it.

Can anyone help me find how can I style this in the Mapbox GL editor (and or via code)?


texts: Using the regexp_substr function in QGIS to get the first letter after white space

texts:
I am looking to get the first letter and second letter (after white space) in the QGIS variable @user_full_name to provide me with the initials of the QGIS user that has created a project.

For instance, James Lewis > JL

My understanding is that the regexp_substr function can be used to achieve this, however I am not sure how the syntax should look to access the first character after the second white space.


texts: KMZ import does not create columns in attribute table

texts:
I would like to import a KMZ file with data on basin internal evaporation recycling into QGIS.

However when I import the KMZ file into QGIS, the relevant data (Bier-jan, BIER-feb, ...) is not added in individual columns of the attribute table but rather all values are added into one single column with together with html syntax.

The cell content of cells in the 'description' column looks is as follows:

<META http-equiv="Content-Type" content="text/html">

<meta http-equiv="content-type" content="text/html; charset=UTF-8">

BIER-jan

0,021422

BIER-feb

0,02088

BIER-mar

0,019949

BIER-apr

0,013789

BIER-may

0,011909

BIER-jun

0,007501

BIER-jul

0,007841

BIER-aug

0,005692

BIER-sep

0,005994

BIER-oct

0,00542

BIER-nov

0,009162

BIER-dec

0,016144

BIER-annual

0,013681

How can I import the KMZ file and add the individual data sets (BIER-jan, BIER-feb, ... ) to dedicated columns in the attribute table?

I already tried to extract the respective numbers with string extraction and string position via the field calculator but a lot cells ended up with a NULL values.
... something like extracting the eight values after the start position of BIER-jan+syntax

 substr(description,strpos(description,'BIER-jan')+17,8)

There is also an online converter which does the job but it is limited to 5MB per month.


texts: Casting geometry to MULTI using GeoPandas?

texts:
Shapefiles can mix simple POLYGONs and MULTIPOLYGONs in the same data source. Spatial databases like PostGIS and SpatiaLite are strict, and will not put a POLYGON in a MULTIPOLYGON geometry column.

I've gotten used to using ST_Multi to fix this issue. But now I am trying to use GeoPandas to do some file processing, including converting from shapefile to GeoPackage (with a bunch of stuff in the middle), and I am running into this error:

gdf.to_file("garbage.gpkg", "GPKG")

ValueError: Record's geometry type does not match collection schema's
   geometry type: 'MultiPolygon' != 'Polygon'

Is there a GeoPandas equivalent to ST_Multi that I can use to fix the geometry before saving to the GeoPackage or SpatiaLite format?


texts: Merging Rasters within Geodatabase using ModelBuilder in ArcGIS Pro?

texts:
I have a geodatabase containing multiple rasters across the same extent with values of either null or 1. The goal is to merge these so where values of 1 overlap, the cell values are aggregated. Exactly the same as '+' in the Raster calculator tool.

This has worked using 'Mosaic to New Raster' however each raster dataset must be input separately into the model. I want to input all the data into the tool from one location.

Using ModelBuilder in ArcPro 2.1.2


texts: Converting BMP to GeoTIFF using GDAL (with aux.xml file)

texts:
I have a BMP with known georeference properties. It is in EPSG:3857, and I have an aux.xml file as well as knowing the lat/lon extent of the image. How do I convert this to GeoTIFF properly?

I have tried the following code which seems to convert the image, but gdal2tiles will not work with it because I think I am missing information.

gdal_translate -of Gtiff -a_ullr -151.8860639 56.77489167 -52.91295 14.56403889 -a_srs EPSG:3857 file.bmp out.tif

Am I missing something there? Below is the aux.xml file. I am not sure exactly how I should use it to properly create the geotif. Is there a way for gdal to open the file and add the projection information? I have also included a public link to Google Drive with the .bmp and aux.xml file, in case that might help.

<PAMDataset>
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</SRS>
  <GeoTransform> -1.6907879268624876e+07,  1.4321638485317408e+03,  0.0000000000000000e+00,  7.7142468591213180e+06,  0.0000000000000000e+00, -1.4321638485317408e+03</GeoTransform>
  <PAMRasterBand band="1">
    <ColorInterp>Gray</ColorInterp>
  </PAMRasterBand>
</PAMDataset>


texts:
I've been on this site blogging about raster files over the week, & this community has proved to be a lot of help

I currently have a map (see below), that I have converted to a raster file, manipulated the map to extract the data I want, and now I just want to save the file, but I keep getting this error message:

"Error in .startAsciiWriting(x, filename, ...) :
  x has unequal horizontal and vertical resolutions. Such data cannot be stored in arc-ascii format"

I have used this code for other images, and it worked just fine.  But for some reason, with this one, I cannot do the last step & save the file.

Does anyone see why this would happen?  And if it would help, I would be happy to post code with an image that did work.

Here is my code:

library(raster)
library(maps)
library(magick)
library(rgdal)

# first need to load working directory

setwd("working here")

# i set this up in a way to use it in a loop later on.

Species_id<-1
Species_i <- image_read(paste("Species (", Species_id, ").jpg", sep=""))

Species_i_rotated<-image_rotate(Species_i,-4) # the -4 is the rotation
image_write(Species_i_rotated, path = "Temporary.png", format = "png")

#Read the temporary image as a raster format using the raster package
#Specify the image band that you are interested in

ConvertedR <- raster("Temporary.png", band = 1)
ConvertedG <- raster("Temporary.png", band = 2)
ConvertedB <- raster("Temporary.png", band = 3)

#isolate color frequency of interest

RGB_selection <- ConvertedR >= 0 & ConvertedR < 50 & ConvertedG >= 0 & ConvertedG < 30 & ConvertedB>=60 & ConvertedB <=90
RGB_selection[100:200,300:400] <- 0 # gets rid of the legend

# Cut the raster to fit as precisely as possible the border of Serbia

ext<-c(51,478,104, 642) # pre-tested in different code.  Double checked and it's fine for what I have
Converted.cr<-crop(RGB_selection, ext)

#Define the coordinates of the borders of the rasters using the sr-org projection 8681 - bosnia.Those need to be adjusted to ensure a final match
# to the data, found @ this website: http://spatialreference.org/ref/sr-org/8681/

xmin(Converted.cr) <- 6588937.517403-20000
xmax(Converted.cr) <- 6907022.282226
ymin(Converted.cr) <- 4668568.705537
ymax(Converted.cr) <- 5117761.487247

#Project the raster into the new coordinate system

crs(Converted.cr) <- "+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel +towgs84=682.0000,-203.0000,480.0000,0.000000,0.000000,0.000000,0.00000000 +units=m +no_defs "
newproj<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
Converted.proj <- projectRaster(Converted.cr, crs=newproj)

#Plot the species map with our adjusted raster
plot(Converted.proj, main="5>25%")
map(region="Serbia", add=T)

# match is great, so save project

writeRaster(Converted.proj,paste("~/Desktop/Heinz_Yan/Final_1",".asc", sep=""), type="ascii", overwrite=TRUE)

# Output
 Error in .startAsciiWriting(x, filename, ...) :
 x has unequal horizontal and vertical resolutions. Such data cannot be stored in arc-ascii format

Here is an example image I am working with


texts: Storing different data types in single column using PostGIS?

texts:
I have three different type of geospatial data Road Direction, Pinpoint location with/without radius, Custom define area which I want to store. If I understood it right this would need three different data types for postgis column, LINESTRING, POINT, POLYGON respectively.

Now the problem for me is these three different types represents a single entity which you can say place_user_interested_in

Can I store these three different types in single column?

The reason why should I store this in single column is to make it more flexible to add more data types later if needed and it would be easier to perform same actions on all of these in single query and some other relevant stuff


texts: d3 + Leaflet apply class value from geojson data

texts:
I have a geojson file that I'm able to draw just fine with D3 + Leaflet. I'm having trouble applying a class to the svg objects with a value from the geojson file. Each feature in the geojson file has a 'handle' attribute. I'm using:

feature.attr("class", function (d) {
    return "code " + d.handle;
});

to try and apply the value, but all of the objects come up with the class of "code undefined". I don't think I'm binding the data correctly..or at all to the svg objects. Any ideas?

bl.ocks.org version:
http://bl.ocks.org/jamierob/a21843bf2b7f3eb54456

js:

<script>
    var map = new L.Map("map", {center: [46.861967, -113.982825], zoom: 17})
            .addLayer(new L.TileLayer("http://tile1.map.umt.edu/tiles/tiles.py/composite/{z}/{x}/{y}.png"));
    var svg = d3.select(map.getPanes().overlayPane).append("svg"),
            g = svg.append("g").attr("class", "leaflet-zoom-hide");
    d3.json("bldgs.geojson", function (collection) {
        var transform = d3.geo.transform({point: projectPoint}),
                path = d3.geo.path().projection(transform);
        var feature = g.selectAll("path")
                .data(collection.features)
                .enter()
                .append("path");
        map.on("viewreset", reset);
        reset();
        // Reposition the SVG to cover the features.
        function reset() {
            var bounds = path.bounds(collection),
                    topLeft = bounds[0],
                    bottomRight = bounds[1];
            svg.attr("width", bottomRight[0] - topLeft[0])
                    .attr("height", bottomRight[1] - topLeft[1])
                    .style("left", topLeft[0] + "px")
                    .style("top", topLeft[1] + "px");
            g.attr("transform", "translate(" + -topLeft[0] + "," + -topLeft[1] + ")");
            //apply a class with the 'handle' value from the geojson file
            feature.attr("class", function (d) {
                return "code " + d.handle;
            });
            feature.attr("d", path);
        }
        // Use Leaflet to implement a D3 geometric transformation.
        function projectPoint(x, y) {
            var point = map.latLngToLayerPoint(new L.LatLng(y, x));
            this.stream.point(point.x, point.y);
        }
    });
</script>


texts: Formula for calculating short distances in UTM?

texts:
I'm fairly new at calculating distances with precision, but I wanted to know what is the best formula for calculating short distances in UTM coordinates.

The distances that I'am calculating are all less than 30 km and all coordinates lie at Central America.

For instance, the pythagorean distance could cause too much bias?

If not, would it still be a respectable choice either way?

If possible, base your answer on some sort manual or academic paper (and make a reference to it).


texts: Why does the 'Raster to vector' conversion give me out lines instead of polygons?

texts:
I've downloaded the Corine Landcover Raster for Europe and wanted to convert the tiff data into polygons.

But when i use the conversion tool, I'll get lines as output, which is useless for me because i need the data of the polygons.

I'm using QGIS 2.10.1.

Here's the data set. I downloaded the 250x250 raster

http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-2


texts:
This is my first time using QGIS.

I downloaded some WFS (http://www.gdos.gov.pl/access-to-geospatial-data) data to use on the Natural Earth QuickStart Kit map (http://www.qgistutorials.com/en/docs/making_a_map.html).

The WFS layers have been added, however they are hidden under the map, rather than being on the map.

How can I fix it?


texts: Can Python 3.6 Script be turned into tool that will work in ArcGIS Desktop?

texts:
I have a Python 3.6 script that I want to turn into a tool for that I can use in ArcMap (10.5). Will that even work or will I strictly have to do it in ArcGIS Pro?

I couldn't seem to find any documentation on this.


texts:
I use PdfFileMerger from PyPDF2 in a plugin that I want to migrate from qgis 2 to 3 but it seems that this lib doesn't exist in qgis 3 (python 3.6).

How I can do this?


texts:
I have an Esri shapefile with following in the .prj file (broken into multiple lines for comprehension):

PROJCS[
  "SGL_Stereo",
  GEOGCS[
    "GCS_SGL_GEO",
    DATUM[
    "D_SGL",
    SPHEROID["Clarke_1880_IGN", 6378249.2,293.4660212936265]],
    PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
    PROJECTION["Double_Stereographic"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",39.15],
    PARAMETER["Scale_Factor",0.9995341],
    PARAMETER["Latitude_Of_Origin",34.2],
    UNIT["Meter",1.0]
]

The problem is that this is a non-standard projection, so any data that I get from the internet, especially that in WGS84, is misaligned, even after re-projecting it.

I found a solution for this, which includes a datum transform. When I try to add a custom datum transform in QGIS (Properties -> CRS -> Datum transformations -> Add), what I see is the following:

I select the source CRS, and match it with WGS84, but I get no further options to enter the parameters of the transformation.

What am I missing?


texts: Export raster layer in python to specific folder

texts:
I am trying to export existing raster layers as tif files. I just want to avoid saving each raster manually. this is the code I am trying to use, but it doesn't seem to work:

from PyQt4.QtCore import *
from qgis.core import *
from qgis.gui import *
import os
import sys

layers = iface.layerTreeView().selectedLayers()
for layer in layers:
    output='Desktop'%layer.name%'.tif'

As specified in a comment below, I modified my code to this:

from PyQt4.QtCore import *
from qgis.core import *
from qgis.gui import *
import os
import sys

layers = iface.layerTreeView().selectedLayers()

for layer in layers:
     file_name = 'D:\tif' + layer.name() + '.tif'
     file_writer = QgsRasterFileWriter(file_name)
     pipe = QgsRasterPipe()
     provider = layer.dataProvider()

     if not pipe.set(provider.clone()):
        print "Cannot set pipe provider"
     continue

     file_writer.writeRaster(
        pipe,
        provider.xSize(),
        provider.ySize(),
        provider.extent(),
        provider.crs())

I can't seem to get it to work and I am not getting any error.

I already selected the layers and it still does not work.


texts: Skip intermediate raster creation in interpolation model

texts:
I am interpolating values at specific coordinates. I am starting with point data and am using ModelBuilder so I can generate interpolated values for each day.

The model (1) iterates through each row (each row is a day), (2) subsets the interpolation to an area within a buffer, (3) creates an IDW raster, (4) samples at the coordinate of interest, and (5) then converts to xls. This works fine, but takes time.

How can I make this more efficient?

Instead of creating a raster and then sampling at a point of interest, how do I interpolate at a specific point?

Thank you.


texts:
I am new to python geospatial programming. i ran the following script and got the corresponding error message

>>> import osgeo
>>> import osgeo.ogr
>>> shapefile = osgeo.ogr.Open("tl_2009_us_state.shp")
>>> numLayers = shapefile.GetLayerCount()

Traceback (most recent call last):
    File "<pyshell#5>", line 1, in <module>
    numLayers = shapefile.GetLayerCount() AttributeError: 'NoneType' object has no attribute 'GetLayerCount'


texts: OpenLayers Map Draw polygon Feature inside boundaries

texts:
Is it possible to draw polygon inside specific boundaries?

draw_interaction = new ol.interaction.Draw({
    source: vector_layer.getSource(),
    type: /** @type {ol.geom.GeometryType} */ ($geom_type.val())
});

map.addInteraction(draw_interaction);

Specify boundaries area in yellow is it possible to add new Geometry Type inside yellow area only? https://jsfiddle.net/vishal05/wb3m90eb/


texts:
I am running on the following setup:

Hardware:

SuperMirco Server
2x Intel Xeon E5-2620 v3 @ 2.40GH
64GB RAM
RAID 5 SATA array with OS, programs and cached tiles
RAID 5 SAS array with Postgres data store

Software:

Windows Server 2012 R2 Standard
Tomcat 8.5.4
PostgreSQL 10.2 with PostGIS 2.4.3
GeoServer 2.9.0 using intergrated GeoWebCache 1.9.0

I have 6 WMS layers which are used by a website (OpenLayers 4.6) and also generate custom images for reports.

For the most part the system is running smoothly but once or twice a week GeoServer stops responding.

No tiles are served from GeoWebCache and direct requests to GeoServer
don't generate an image.
PostgreSQL is still running with no issues.
I can log into the GeoServer control panel.
I have tried restarting GeoServer from the Tomcat Web Application Manager and this doesn't solve the issue.
The only solution I have found is to restart Tomcat.

Any suggestions and/or where to start looking for the problem?