texts: ArcObjects ExportMapImage modifies image extent

texts:
I am trying to extract a map image at a certain range. I am using IMxdServer's ExportMapImage method, documented here.

Here is the relevant snippet of code:

    IMapDocumentPtr ipMapDoc(CLSID_MapDocument);
    IMapDescriptorPtr ipMapDescriptor(CLSID_MapDescriptor);
    IMxdServerPtr ipMapServer(CLSID_MxdServer);
    IArrayPtr ipArray(CLSID_Array);
    IUnknownPtr ipUnknown;
    int dpi = 300;
            IEnvelopePtr ipPixelBoundsEnv(CLSID_Envelope), ipReturnExtent(CLSID_Envelope);
    ipPixelBoundsEnv->PutCoords(mUpperLeftLong, mLowerRightLat, mLowerRightLong,  mUpperLeftLat);

    ipMapServer->Start(MxdPath);
    ipMapServer->get_MapDescriptors(&ipArray);

    ipArray->get_Element(0, &ipUnknown);
    ipMapDescriptor = (IMapDescriptorPtr) ipUnknown;

    IColorPtr ipColor(CLSID_RgbColor);
    ipMapDescriptor->get_BackgroundColor(&ipColor);
    ipMapDescriptor->put_TransparentColor(ipColor);
    ipMapDescriptor->put_Extent(ipPixelBoundsEnv);
    ipMapDescriptor->put_FullExtent(ipPixelBoundsEnv);
    ipMapDescriptor->put_KeepAspectRatio(VARIANT_FALSE);

    HRESULT hr = ipMapServer->ExportMapImage(ipMapDescriptor, dpi, (256*(hortizonalTileBlockCount+2)), (256*(verticalTileBlockCount+2)), CComBSTR("PNG8"), bstrFileLocation, &ipReturnExtent);
    double ux, uy, lx, ly;
    ipReturnExtent->get_XMax(&ux);
    ipReturnExtent->get_XMin(&lx);
    ipReturnExtent->get_YMax(&uy);
    ipReturnExtent->get_YMin(&ly);
    std::cout<<"return extent "<<lx<<","<<ly<<","<<ux<<","<<uy<< " \n";

In my current example the following bounds are passed in:

-75.9375,39.6395,-77.3438,40.714

The return extent is shifted a touch:

-77.3438,39.4736,-75.9375,40.8799

I'm extracting quite a few images but I put one example extent for reference.


texts:
I have two layers, one with lines, one with polygons. I want to delete lines that fully fall (contain) within the polygons. It seems really simple, but I can't figure this one out in either PostgreSQL/PostGIS or QGIS.
For the SQL I am trying:
delete from lines_test2
 where ( select a.geom, b.geom
 from lines_test2 a, straten b
 where st_intersects(a.geom, b.geom));

But it gives me the error: 'subquery must return only one column'. Clipping in QGIS isn't a viable option either, because it should only delete full lines not parts of them.


texts: Geocoding adresses under privacy concerns

texts:
I have a database of home addresses and service locations. People drive to the service locations. I'd like to create a map with transparent colored circles indicating the average drive times people face getting to our service locations.

I know I can easily get the drive times using Google's Developer API, but I have a slight privacy concern: I'm not sure we're allowed to use the addresses in a way that would necessarily share them with google.

So, the question is really this: is there a way I can geocode these locations and then calculate drive times between them, without having to use a service that will store the data?

EDIT: Sorry for the lack of clarity. This is in the US and I have two basic tasks: (1) geocode addresses, and (2) get driving times between two addresses. In each case, I am looking for a software/tool that does not require disclosing the addresses to an outside source. My understanding is that we can use the datasciencetoolkit to geocode the addresses locally. But, I don't know what I can use to get directions.


texts:
I have a Enterprise Geodatabase (SDE) in SQL Server 2016.
I have a nightly python script to preform some maintenance including the compress function (among other things).

I am noticing that the state counts after the compress are not going lower than 8. The before state count is higher but no matter what the count never reaches 1 (sometimes it does go to 7).

Why is the count staying at this number after the compress? does this indicate something is wrong, or a possible conflict/issue in the near future?


texts:
In our application, we exposed a custom file format to ArcScene through a custom plugin data-source.
However, we experience some problems which we ascribe to caching issues.

We generate a dataset with 5 fields, in particular OID, Shape and 3 custom columns. This dataset is added
to the scene through a standard feature layer. ArcScene properly calls IPlugInDatasetHelper.get_Fields to obtain the IFields from our dataset.

Upon user interaction, this dataset might get regenerated where the number of fields might change.
In this case to OID, Shape and a single custom row.
If ArcScene now calls IPlugInDatasetHelper.FetchAll ,
the last parameter (object FieldMap) still contains 5 entries. ArcScene does not call IPlugInDatasetHelper.get_Fields anymore, despite the column count having changed.

Is there a way to tell ArcScene/ArcMap the fields in our custom data-source have changed?

We already tried ISceneGraph.Invalidate on our feature layer or ActiveView.Redraw of no avail.

Thanks!


texts: Is it possible to convert between shapefile and geojson using geotools without access to file system? I have that restriction in my web app, but I need to create shapefile as byte array and then do something with that. Unfortunately I only found a way where I need to specify new data store based on file path.


texts: Autoincrement field when New Feature added

texts:
I know that this is probably a very basic programming but I am brand new to Python.

I have a layer in a geodatabase that I need a unique ID field (for joining). I can't use the ObjectID because that can change. I can't use a GlobalID because the data will be used by databases that can't handle the "{" brackets or dashes.

I have already populated the field. When I add a new feature I want to have the value of the max existing value + 1 for the new feature.

Thanks to Jason here is what I have so far-

import arcpy
gdb_full_path = r'D:\GIS\ZA.gdb\Current'
feat = 'ZA_Merge'
arcpy.env.workspace = gdb_full_path
with arcpy.da.SearchCursor(feat,["ANID"]) as cursor:
     for row in cursor:
         unique_id_counter = row.count
with arcpy.da.InsertCursor(feat,["ANID"]) as cursor:
     cursor.insertRow((unique_id_counter + 1))

Here is the error that I get -

Runtime error
Traceback (most recent call last):
  File "", line 9, in
TypeError: unsupported operand type(s) for +: 'builtin_function_or_method' and 'int'


texts: Transform / rectify unprojected raster using target extent co-ordinates

texts:
I have many unreferenced raster tiles that I am trying to procedurally georectify. For each tile I know the exact extent co-ordinates (xmin, xmax, ymin, ymax) in my target CRS.

How can I transform these rasters, translating their arbitrary pixel-based co-ordinates (e.g 0,0 | 0,1280 | etc.) to match those calculated from my extent values?

All the GRASS tools seem to want manually generated GCPs and I'm unsure which GDAL tools would be suitable, given the raster tiles are unprojected.

I am looking for python / command line solutions on OSX.


texts: raster::plot() plotting raster too small- is there an alternative on ggplot?

texts:
I am trying to plot this raster on R. The object is quite big, but when plotting it using the raster::plot() function, I obtain a very tiny image, as attached.

The code I use to produce it is the following:

r <- raster::raster("MOD13Q1_2000-02-18.250m_16_days_NDVI.tif")
r <- setMinMax(r)
raster::plot(r)

I tried to set xlim limits, with no success. I then tried to move to ggplot2, as an alternative (following here), but I obtain the following

ggplot(r) + geom_tile(aes(fill=value))
Error: ggplot2 doesn't know how to deal with data of class RasterLayer

I found around that, to plot raster objects with ggplot(), I should first convert it into a data.frame. However, as my object is quite big (about 55million obs.), that would be very inefficient. I am then stuck on what to do. Any thoughts?


texts:
ArcMap 9.3.1 w/ ArcFM.  Using C#.

The goal is to be able to prevent certain edits based on specified business rules.

One option might seem to be looking at all operation edits in BeforeStopEditing and somehow cancelling the operation here.  This is not possible as ESRI states you can't call AbortOperation here and throwing an exception does not prevent the operation from committing.  Even setting a flag in BeforeStopEditing and calling UndoOperation in StopEditing if the flag is set won't work.

The obvious documented approach is to create a custom class extension and apply to all (required) tables.  A given table can only have one class extension applied.  Since we are using ArcFM it has applied its class extension(s) to the feature class tables.  We could probably use this approach for the data tables.

We are currently using the Editor IEditEvents* to listen for feature edits, but this doesn't provide data table edits.

We are also hoping to avoid applying a custom class extension to 100's of tables, the maintenence group would not agree to that.

It's also not feasible to get all IObjectClasses and listen to each directly in the code.

The hope was to find an existing coclass or interface that broadcasts edit events for all tables.  The doc's seem to suggest ObjectClassEvents may do this, but only in VB.  Trying to create an instance of this in C# throws an exception (factory can't create).

Suggestions ?


texts: Convert raster into array within a PostGIS database

texts:
I need to convert a raster (e.g. GeoTIFF format) into an array (it can be any kind of array e.g. numpy array).
So far I managed to do that in a python script (as they describe it at this page https://gist.github.com/gzagatti/a649674d85b480b98b7d75b31d4c9555).
Now I am trying to do this within the PostGIS database. Is there a method or function from PostGIS which can convert the raster directly into an array?
I'm using PostgreSQL 12 with the PostGIS 3.


texts: Displaying data from standalone PostGIS using ArcGIS API for JavaScript?

texts:
I'd like to serve out data from a PostGIS server and consume it as a layer in the ArcGIS JavaScript API. Searching for "arcgis javascript postgis" I see information on ArcGIS Desktop, ArcGIS Server and information on creating a FeatureLayer.

My experience with Esri is with ArcSDE from about 10 years ago. At this stage, I'm left with no obvious clues. Browsing the documentation, I found this page on Layers and Data. Looking at the Core Layer Types section, it looks like a FeatureLayer isn't what I'm looking for - my database will be on Amazon RDS, not in the ArcGIS cloud.

Is there a way to configure a layer in JavaScript that will point to my PostGIS server or is it a little more involved than that?

For example, I'm thinking that I could create a web service call that returns the data from the PostgreSQL database. I could then consume this data and generate a GraphicsLayer.


texts:
I have a fairly large geotiff that has only a small area of interest. I have used clip by rectangle to generate a new geotiff of this area of interest but the resulting geotiff is the same size as the original raster (in both size and extent), with the area of interest present but surrounded by nodata.

Is there way to essentially crop the rectangle of real data out from the nodata?


texts: How to run an Identify task on a WMS layer

texts:
I'm new to WMS layers, and am looking for a way to retrieve a layer's attributes when clicking on the map, using the ArcGIS Server JS API.

The documentation is a bit sparse, but from this question it appears that I'll need to construct a query using GetFeatureInfo.

Before I re-invent the wheel, has anyone written a wrapper function or sample code which builds the appropriate query when the user clicks on the map?


texts:
I have a software package that requires "Unprojected GeoTIFFs" and no matter what I try I cannot get them accepted. Ive been trying GDAL but to no avail. I have an example of a working file and one I have exported ie that I want to use so can see th edifference, just dont know how to modify it, any advice? Here is the gdalinfo outputs,

  i:\PNG_DATA>gdalinfo workingexample.tif Driver: GTiff/GeoTIFF Files:
  workingexample.tif Size is 1624, 1083 Coordinate System is:
  GEOGCS["WGS 84",
      DATUM["WGS_1984",
          SPHEROID["WGS 84",6378137,298.257223563,
              AUTHORITY["EPSG","7030"]],
          AUTHORITY["EPSG","6326"]],
      PRIMEM["Greenwich",0],
      UNIT["degree",0.0174532925199433],
      AUTHORITY["EPSG","4326"]] Origin = (10.477938743014445,59.473562724738137) Pixel Size =
  (0.000007736964863,-0.000003930129367) Metadata:   AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)   TIFFTAG_XRESOLUTION=96
  TIFFTAG_YRESOLUTION=96 Image Structure Metadata:   INTERLEAVE=PIXEL
  Corner Coordinates: Upper Left  (  10.4779387,  59.4735627) (
  10d28'40.58"E, 59d28'24.83"N) Lower Left  (  10.4779387,  59.4693064)
  ( 10d28'40.58"E, 59d28' 9.50"N) Upper Right (  10.4905036,
  59.4735627) ( 10d29'25.81"E, 59d28'24.83"N) Lower Right (  10.4905036,  59.4693064) ( 10d29'25.81"E, 59d28' 9.50"N) Center      (  10.4842212,  59.4714346) ( 10d29' 3.20"E, 59d28'17.16"N) Band 1 Block=1624x1 Type=Byte, ColorInterp=Red Band 2 Block=1624x1 Type=Byte,
  ColorInterp=Green Band 3 Block=1624x1 Type=Byte, ColorInterp=Blue Band
  4 Block=1624x1 Type=Byte, ColorInterp=Undefined

  i:\PNG_DATA>gdalinfo Overviewtiff.tif Driver: GTiff/GeoTIFF Files:
  Overviewtiff.tif Size is 1536, 928 Coordinate System is: PROJCS["WGS
  84 / UTM zone 55S",
      GEOGCS["WGS 84",
          DATUM["WGS_1984",
              SPHEROID["WGS 84",6378137,298.257223563,
                  AUTHORITY["EPSG","7030"]],
              AUTHORITY["EPSG","6326"]],
          PRIMEM["Greenwich",0],
          UNIT["degree",0.0174532925199433],
          AUTHORITY["EPSG","4326"]],
      PROJECTION["Transverse_Mercator"],
      PARAMETER["latitude_of_origin",0],
      PARAMETER["central_meridian",147],
      PARAMETER["scale_factor",0.9996],
      PARAMETER["false_easting",500000],
      PARAMETER["false_northing",10000000],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
      AUTHORITY["EPSG","32755"]] Origin = (500200.000000000000000,9258575.000000000000000) Pixel Size =
  (25.000000000000000,-25.000000000000000) Image Structure Metadata:
  COMPRESSION=LZW   INTERLEAVE=PIXEL Corner Coordinates: Upper Left  (
  500200.000, 9258575.000) (147d 0' 6.51"E,  6d42'27.35"S) Lower Left  (  500200.000, 9235375.000) (147d 0' 6.52"E,  6d55' 2.88"S) Upper Right (  538600.000, 9258575.000) (147d20'57.34"E,  6d42'26.90"S) Lower Right (  538600.000, 9235375.000) (147d20'57.89"E,  6d55' 2.41"S) Center      (  519400.000, 9246975.000) (147d10'32.07"E,  6d48'45.00"S) Band 1 Block=1536x1 Type=Byte, ColorInterp=Red   Mask Flags: PER_DATASET
  ALPHA Band 2 Block=1536x1 Type=Byte, ColorInterp=Green   Mask Flags:
  PER_DATASET ALPHA Band 3 Block=1536x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA Band 4 Block=1536x1 Type=Byte,
  ColorInterp=Alpha

  i:\PNG_DATA>

Anyone got any idea on what command to use to remove the projection info?


texts: Merging fields when doing merge of different layers

texts:
I'm trying to merge 6 points layers, each containing the same fields (Scientific Name, Lat, Long, Count, Common Name) but for different provinces. When I merge the layers, in the geoprocessing wizard, the fields come up in field map as separate for each layer, even though they have the same name. When I completed one of the merges, there would be 2 layers for Common name, for example, with some of the records having "null" in one of the fields and then having data in the other Common Name field. For count, however, this doesn't seem to be a problem.

How do I merge the layers and fields so that all the data from the layers is in just 1 field instead of multiple fields?


texts: customer prospecting for the united states within specific drive time

texts: I have a set of demographic variables that I want to use to identify areas within United States that meet this criteria. I have used the customer prospecting tool in Business Analyst but the issue that I have is that I want to identify those areas only within 10 min of a drive time.


texts: How to focus the image when it is added to the map in OpenLayers?

texts:
I created a map with OpenLayers OSM. On this map I add an image provided by GeoServer. My question is how to focus the image when it is added to the map. As zooming into the area on the map where the image is added.

var result = new ol.layer.Image({
    source: new ol.source.ImageWMS({
        url: 'http://localhost:8080/geoserver/pgrouting/wms',
        params: params
    })
});
map.addLayer(result);

map.zoomToMaxExtent();


texts: Distance difference when moving longitude in maps services

texts:
I try to get two juxtaposed tiles in maps services like Google Maps or Bing Maps at a given zoom level (20 in my case). I managed to do it with those methods:

def meters_per_pixel_NS(lat, zoom):
    return 156543.03392 * math.cos(math.radians(lat)) / 2**zoom

def meters_per_pixel_EW(zoom):
    return 156543.03392 / 2**zoom

def get_next_coord(coord, dist):
    K = 1852*60 # 1 nautical mile  = 1852 m; multiplied by 60 to get degrees
    delta = dist/K
    return coord + delta

Let's say I have the following tile of 256 pixels, captured in 59.329411, 18.059484 somewhere in Stockolm and the following tile that would be juxtaposed on the right, captured in 59.329411, 18.05982793911214

According to my method meters_per_pixel_EW(zoom), the distance from the left border to the right border (which is also the shifting) is 38.22 meters. It is a constant that does not depend on the latitude and which is wrong, because it is actually 19.57 meters wide, just like the value I would get if I did a vertical shifting (that depends on the latitude).

That seems logical, because otherwise the image would be completely distored (a height of 19.57 meters and a width of 38.22 meters).

And yet, I managed to obtain the correct coordinates for the next tile by converting 38.22 meters in degrees. Again, it seems logical as the mercator projection does not distort the longitudes.

In short, I asked a shifting of 38.22 meters and I received a shifting of 19.57 meters. If I did the same in Mexico, I would ask the exact same shifting (in degrees) and get a higher horizontal distance.

I checked both on Google Maps and Bing Maps, it works. How can I explain this "paradox"? Is there something I don't understand in Mercator projection or is there a transformations/corrections that is made by maps services ?


texts: Dealing with Very Large Satellite Imagery

texts:
I'm curious if anyone has suggestions for viewing and processing really large imagery. I'm pretty experienced with desktop apps like ENVI and am comfortable with python and GDAL but I'm analyzing dozens and dozens (soon to be hundreds) of DigitalGlobe images that are around 40 gigs per file and I'm struggling to effectively do data exploration because it takes so long to even view them. And then I'm not sure what the best environment is to handle hundreds of images this size.

Does anyone know any cloud based viewers and analysis sandboxes that you recommend? I'm getting access to some beefy machines over the next week or two and I'd prefer to do it locally. Overall I just wanted to get an idea of suggested pipelines for analyzing this much imagery.


texts: I am creating a sea chart. The sea chart has ships lanes in a line layer. In the same layer, the safe depth of the ship lane is found in a column. I want to use this safe depth as a label on the ship lane, which works without a problem. But to follow standard themathising rules and get a clear view of the number, I need the ships lane line not to be visible underneath and just around the number (see as an example of what I want QGIS to do the number 2,4 on the image below). Creating a white background under the text is not an option, this could cover other important features.  So is it possible to make the line just disappear under the label and for instance 2mm on the sides of the label?


texts:
I am trying to migrate from serving static MVT to generated on the backend by PostGIS with ST_AsMVT(). Or find a way of minimizing the time of building MVT from the database.

I have a dataset of political, historical maps, and I am serving them as prerendered Mapbox Vector Tiles. Dataset itself is not massive. It's 2300 entities for 200 years ~ 1Gb of data. As a result, distant zoom tile (/1/0/0) has ~ 600 geometries. They are filtered on client side with filter option by year in MVT style. We want to introduce an ability for the end-user to upload their maps.

Pipeline to produce an MVT tileset that we are using right now:

Loop over geometry in the database and dump it in a file as a GeoJSON feature (feature collection does not fit into memory)
Process output with tippecanoe and zoom up to 10
Serve file.mbtiles with mbtiles-server

It took around 30 minutes and full CPU load for the whole cycle.

As one way of optimization here - we could rebuild only the part that was changed and update tileset with tile-join

I've tried to serve MVT directly from PostGIS with ST_AsMVT(),
Geometry is in SRID:4326 queries were executed with similar simplification level and against the same bounding box /1/0/0

Full queries and explanation are here

The straight-forward approach of MVT generation took 130-177 seconds
Simplification after transforming to SRID:3857 (WebMercator) took 30 seconds
Simplification before processing to SRID:3857 took 3 seconds

-- three seconds to execute
SELECT ST_AsMVT(a, 'stv')
FROM (
  SELECT
      api_spacetimevolume.id
      , api_spacetimevolume.start_date::INTEGER
      , api_spacetimevolume.end_date::INTEGER
      , api_spacetimevolume.references
      , ST_AsMVTGeom(
          ST_SnapToGrid(
            ST_Transform(
              ST_Simplify(
                  api_spacetimevolume.territory
                , 0.44 -- Dynamic based on zoom level
              )
              , 3857
            )
            , 1
          )
          , TileBBox(1, 0, 0)
      ) as territory
      , api_spacetimevolume.entity_id
      , api_territorialentity.wikidata_id
      , api_territorialentity.color
      , api_territorialentity.admin_level
  FROM api_spacetimevolume
  JOIN api_territorialentity
  ON api_spacetimevolume.entity_id = api_territorialentity.id
  WHERE territory && TileBBox(1, 0, 0, 4326)
) as a
;

Even with 3 seconds per tile, it is obvious that it should be cached on disk or in the database. Temporary caching (nginx or redis) will not solve the problem completely.

In this article Michal Zimmerman is generating tiles in advance by PostGis and storing them on disk.

Well, the questions

What do you recommend to speed up the tile delivery?
What is the best way to keep the tile generation flexible?


texts: Map projections for each quadrant of the globe

texts:
I am using an xarray data and matplotlib to plot it before export
This is the xarr data used
<xarray.DataArray 't2m' (latitude: 361, longitude: 720)>
array([[273.9, 273.9, 273.9, ..., 273.9, 273.9, 273.9],
       [274.1, 274.1, 274.1, ..., 274.1, 274.1, 274.1],
       [274.5, 274.5, 274.5, ..., 274.5, 274.5, 274.5],
       ...,
       [220.6, 220.6, 220.6, ..., 220.6, 220.6, 220.6],
       [221.6, 221.6, 221.6, ..., 221.6, 221.6, 221.6],
       [223. , 223. , 223. , ..., 223. , 223. , 223. ]], dtype=float32)
Coordinates:
  * latitude    (latitude) float64 90.0 89.5 89.0 88.5 ... -89.0 -89.5 -90.0
  * longitude   (longitude) float64 -180.0 -179.5 -179.0 ... 178.5 179.0 179.5
    step        timedelta64[ns] 03:00:00
    valid_time  datetime64[ns] 2020-09-14T03:00:00
Attributes:
    standard_name:  air_temperature

The longitude range is from -180 to +180. If I use the below data with projections the image produced is correct.
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

imgDpi = 100
fig = plt.figure(figsize=(512/imgDpi, 512/imgDpi), dpi=imgDpi) #to convert from inches to px
ax = fig.add_subplot(1,1,1, projection=ccrs.epsg(3857))
ax.imshow(xarr, transform=ccrs.PlateCarree(), origin='upper', cmap='viridis')
ax.axis('off')
fig.tight_layout(pad=0)
plt.show()
fig.savefig('unsliced.png', format="png", transparent=True, dpi=imgDpi)

But if I slice the data into four quadrants and plot, there is an issue with projection
topLeft = xarr.sel(longitude=slice(-180, 0), latitude=slice(90, 0))
topRight = xarr.sel(longitude=slice(0, 180), latitude=slice(90, 0))
bottomLeft = xarr.sel(longitude=slice(-180, 0), latitude=slice(0, -90))
bottomRight = xarr.sel(longitude=slice(0, 180), latitude=slice(0, -90))

axTL = fig.add_subplot(2,2,1, projection=ccrs.epsg(3857))
axTL.imshow(topLeft, transform=ccrs.PlateCarree(), origin='upper', cmap='viridis')
axTL.axis('off')

axTR = fig.add_subplot(2,2,2, projection=ccrs.epsg(3857))
axTR.imshow(topRight, transform=ccrs.PlateCarree(), origin='upper', cmap='viridis')
axTR.axis('off')

axBL = fig.add_subplot(2,2,3, projection=ccrs.epsg(3857))
axBL.imshow(bottomLeft, transform=ccrs.PlateCarree(), origin='upper', cmap='viridis')
axBL.axis('off')

axBR = fig.add_subplot(2,2,4, projection=ccrs.epsg(3857))
axBR.imshow(bottomRight, transform=ccrs.PlateCarree(), origin='upper', cmap='viridis')
axBR.axis('off')


texts: Cannot load layers into ArcGIS Viewer for Flex

texts:
I am trying to load some files into ArcGIS Viewer for Flex. My files are loaded onto the server correctly but trying to publish it into the Flex Viewer is unsuccessful.

First when I search for my layers from the server I get this warning that the program couldn't find crossdomain.xml. (see image below)

I continue to select my layers but when it tries to load I have been getting this error.

UAE Web Map layer failed to load: Fault code: Channel.Security.Error

Fault info: Security error accessing url

Fault details: Destination: DefaultHTTP


texts: Convert lineString to MultilineString in PostGIS

texts:
How can I convert the line string to multiline string in PostGIS?
I have converted geometry points to line string by this: st_makeline(location.shape)
but by this command:

  SELECT uuid_generate_v4()                                         AS uuid,
          j."deviceId",
         st_astext(st_linemerge(st_union(j.shape))) AS multiLine
  FROM j
  group by date, j."deviceId" ;

I got this results:

It's result is again line string?


texts:
By using the following code I am trying to plot Landsat 8 RGB image onto the folium maps. Following is the code:
import rasterio as rio
import folium
from pyproj import Transformer

## LC08 RGB Image
in_path = 'RGB.TIF'

dst_crs = 'EPSG:4326'

with rio.open(in_path) as src:

    img = src.read()

    src_crs = src.crs['init'].upper()
    min_lon, min_lat, max_lon, max_lat = src.bounds

## Conversion from UTM to WGS84 CRS
bounds_orig = [[min_lat, min_lon], [max_lat, max_lon]]

bounds_fin = []

for item in bounds_orig:
    #converting to lat/lon
    lat = item[0]
    lon = item[1]

    proj = Transformer.from_crs(int(src_crs.split(":")[1]), int(dst_crs.split(":")[1]), always_xy=True)

    lon_n, lat_n = proj.transform(lon, lat)

    bounds_fin.append([lat_n, lon_n])

# Finding the centre latitude & longitude
centre_lon = bounds_fin[0][1] + (bounds_fin[1][1] - bounds_fin[0][1])/2
centre_lat = bounds_fin[0][0] + (bounds_fin[1][0] - bounds_fin[0][0])/2

m = folium.Map(location=[centre_lat, centre_lon],
                   tiles='Stamen Terrain', zoom_start = 10)

# Overlay raster (RGB) called img using add_child() function (opacity and bounding box set)
m.add_child(folium.raster_layers.ImageOverlay(img.transpose(1, 2, 0), opacity=.7,
                                 bounds = bounds_fin))

# Display map
m

After running the above code I am getting the following output:

Can someone please let me know how to plot the raster after removing the no-data values i.e. LC08 RGB image without black surrounding borders.


texts:
I am modifying a OpenLayers setup which uses WFS capabilities to create a vectorLayer on a map to support multipoints.

GeoServer is setup to send points and multipoints however only points are displayed. Is there a specific setting that prevents multipoints from being displayed in OpenLayers?

MySet.config.myLayer = {
styleMap: new OpenLayers.StyleMap(
    {
        'default': {
            cssClass: "${categorie}",
            label: "${count}",
            popupActive: "${active}"
        }
    }
),
protocol: new OpenLayers.Protocol.WFS({
    url: MySet.config.WFSurl,
    featureType: "mypoints",
    geometryName: "location",
    version: "1.1.0",
    method: "GET",
    outputFormat: "GML2"

}),
maxResolution: 6.720,
minResolution: 0.210,
renderers: [MySet.Renderer.Anchor],
rendererOptions: {iconSize: MySet.config.iconSize, imgPath: MySet.config.imgPath, classificationInfo: MySet.config.classificationInfo},
strategies: [
    new MySet.Strategy.ThemeVisibility({attributeName: 'categorie'}),
    new MySet.Strategy.Sort({sortHeight: MySet.config.sortHeight}),
    new MySet.Strategy.Cluster({attributeName: 'categorie', attributeValue: 'cluster', threshold: 2, distance: 21}),
    new OpenLayers.Strategy.BBOX()
]};

Draw

createVectorLayers : function () {
            vectorLayer = new OpenLayers.Layer.Vector(null, this.config.myLayer);
            map.addLayer(vectorLayer);
            var keys = ['prov', 'gem', 'wk'];
            for (var i = 0; i < keys.length; i++) {
                var options = this.config.clusterLayer[keys[i]];
                var lyr = new OpenLayers.Layer.Vector(options.title, options);
                map.addLayer(lyr);
            }
        }


texts: Unable to connect to arcgis server local

texts:
I have ArcGIS Server 10.1 installed and running , and published the services its worked , and setup & Configuration Web Adaptor

when i go the ArcGIS Web Applications for the Microsoft .NET Framework to create new web

i cannot add arcgis server local when i type the host machine Name or "localhost" this message showing

Error: Unable to connect to server

what i can do please


texts:
I'm not used to QGIS, and I'm working on a project that involves hundreds of .tiff rasters (interpolations) of attributes of the same place. (the rasters are not really related by time)
In order to export one single raster as a complete map I should click for 300 times on: Project > New print Layout > Put a map, a legend, a title and an arrow > Export as images.
Is there a way to automatically export every raster with his own legend and title (obtained from the name of the associated raster)?


texts: How to calculate purchasing power on postcode level?

texts: Does anyone know of an algorithm to calculate purchasing power on a postcode level? What source data is needed? Purchasing power on admin. level, population?


texts: Change units from degrees to meters in QGIS 3.6

texts:
I'm cleaning up shapefiles (mainly for US areas) in QGIS 3.6 to import into other CAD environments such as Rhino. One problem encountered is that the units are in degrees, resulting in imports being several orders of magnitude wrong. How do I change the unit to meters?

I searched and understand it's something to with the CRS (switching to UTM instead of WGS), but as there seems to be many options/zones, I'm not really sure which is the correct one.


texts: Which software converts TIFF to TJP (Tiled Jpeg)?

texts:
I have Tiff images and I need to export it to Tiled Jpeg. I don't know which software converts TIFF to TJP. Can anyone kindly tell me the name of the software which converts?

Thanks


texts: GeoServer 2.4.1: How to use JSONP in a WFS request?

texts:
I've ready that I can enable JSONP by setting ENABLE_JSONP in web.xml:

  <context-param>
    <param-name>ENABLE_JSONP</param-name>
    <param-value>true</param-value>
  </context-param>

And I've read you can set outputFormat=json&format_options=callback:myCallback but I only get json back.  Is this functionality broken/reverted?  I'm on Windows 8.1.


texts: Problem installing QGIS Valmeira on Windows 7 Professional 64-bit

texts:
I have been having trouble installing QGIS Valmeira on Windows 7 Professional with SP1. When I installed to the default directory C:\Program Files\QGIS Valmiera\, many of the predefined coordinate projections were not loaded. Several other projections were listed in French.

At Andre's suggestion I then reinstalled to a folder without blanks such as C:\QGIS_Valmeira. This worked fine and brought in the previously missing projections.

I did not have this problem with earlier versions of QGIS like Lisboa. What is going on here? Is this a bug in version 2? If so, it must be causing problems for many other Windows users of QGIS.


texts: How to use a secondary geometry in LABEL object

texts:
I have a set of settlement features with two geometries: one is an outline polygon and another is a label point.

The features are fetched from a PostGIS db:

LAYER
    NAME "settlements"
    CONNECTIONTYPE postgis
    CONNECTION "user=x password=x dbname=x host=localhost"
    DATA "polygon_geom from (select id, name, polygon_geom, label_geom from settlements) as subquery using unique id"

The polygon_geom is then used in all CLASS STYLE and LABEL objects of this layer.

Is it possible to make the LABEL objects use the secondary label_geom geometry?

The only solution I've found was to setup a separate layer which gets that label_geom as the primary geometry. But is it possible to achieve it within one layer?


texts: How to create polygons with a set area in MapInfo?

texts:
I need to create polygons with set areas in MapInfo. There is a solution for QGIS in the form of the the QAD plugin but I can't find an equivalent tool for MapInfo.

The task I'm actually trying to do is create polygons which are e.g. 10% the size of an existing polygon.


texts: Choropleth Map: data source by location name (without geometries)

texts:
I am interested in creating a choropleth map using Mapbox, while the data source I have in hand is built only from the location name and the "intensity" value. i.e.:
[{"country": "fr", "city": "paris", "value": 7000}, ... ]

All the code examples and explanations I could get require you to have geographies/coordinates of polygon.
However, I don't have this data, or actually, I have too many cities that getting the polygon for each is going to increase the size of my data source.
Can I actually project data on such map with only country/city name in my data source?


texts:
I am looking to create a predictive model in ArcGIS.  I have done this previously with glm with poisson distribution and found it relatively straightforward in raster calculator e.g.

exp(intercept + (xvariable) + (yvariable))

However, it is likely that i shall need to use a more complex approach to account for both local and landscape variables and will need to use a Generalised Linear Mixed Model.

Does anyone have experience of this?  Is it possible to do in ArcGIS(raster calculator), how would the equation likely differ to that above?


texts: Box cropping shapefile based on particular latitude and longitude using QGIS?

texts:
I have a shapefile and I need to crop/clip it based on a rectangle of particular latitude and longitude using QGIS.

If it is possible, how do I do that?


texts: Filtering interactive leaflet-map with dropdown-menu?

texts:
I try to set up a filterable "shot chart" like http://graphics.latimes.com/kobe-every-shot-ever/ with leaflet. I've a geojson file thats stores all the shots data, I've them showing up in a layer.

GEOJSON-Example:

{
  "type": "Feature",
  "properties": {
    "ID": 2,
    "SEASON_SHO": "1998-1999",
    "GAME_TYPE": "Regular",
    "PERIOD": 1,
    "MINUTES_RE": "06:25:00",
    "ACTION_TYP": "Jump Shot",
    "SHOT_TYPE": "2-PT-Wurf",
    "SHOT_DISTA": 15,
    "GAME_DATE": "05.02.1999",
    "HEIM": "SEA",
    "GAST": "DAL",
    "GEGNER": "Seattle SuperSonics",
    "SHOT_ATTEM": 1,
    "SHOT_MADE_": 0
  },
  "geometry": {
    "type": "Point",
    "coordinates": [
      -0.314959,
      -0.159125
    ]
  }
}

I've set up an layer for all shots and a layergroup to filter the markers.

var url = 'data/Wuerfe.geojson';
var allewuerfe;

allewuerfe = L.geoJson(null, {

    onEachFeature: onEachFeature,
    style: style,
    pointToLayer: function (feature, latlng) {
        return L.circleMarker(latlng, allewuerfe);
    }
});

$.getJSON(url, function (data) {
    allewuerfe.addData(data);
});

var data = L.layerGroup([
    ]);
data.addLayer(allewuerfe);
data.addTo(map);

I already have some filters-buttons for single events, like "reset".

<button class="btn reset-btn" id="all">Reset</button>

document.getElementById("all").addEventListener('click', function (event) {
    data.clearLayers();
    data.addLayer(allewuerfe);
    data.addTo(map);
});

But what I'd like to set up is an dropdown-filter for "shots missed", "shots hit" and "all shots".

           <select name="shot-result" class="custom-select" id="shotresult">
                <option value="all">All Shots</option>
                <option value="Made">Made</option>
                <option value="Missed">Missed</option>
            </select>

The data-info for that is stored in the geojson at "SHOT_MADE_", 1 means shot was made, 0 shot was missed.

I would like to get an code-example of how to get the filter work.


texts:
I have about 600 folders with shapefiles always having the same filename per folder (which contain different data, of course). I've created one QGIS project file to open those shapefiles per folder always in the same way. I did this on one set of shapefiles, adjusted everything to my liking in QGIS and saved it. Afterwards in the code, I've changed the absolute path to those shapefiles to a relative path, and copied in every folder such a qgs-file.

Is it possible to make each project zoom to its extent on opening? Now, every project file I open zooms to the location of my initial "project", where I created the first project file.

I see this in the project file:

    <mapcanvas annotationsVisible="1" name="theMapCanvas">
    <units>meters</units>
    <extent>
      <xmin>57917.75605197192635387</xmin>
      <ymin>175589.96959768730448559</ymin>
      <xmax>66470.8240169670607429</xmax>
      <ymax>188067.17660002846969292</ymax>
    </extent>

or this:

    <projectlayers>
    <maplayer autoRefreshTime="0" refreshOnNotifyMessage="" geometry="Line" readOnly="0" simplifyDrawingHints="1" autoRefreshEnabled="0" refreshOnNotifyEnabled="0" simplifyAlgorithm="0" labelsEnabled="0" minScale="1e+8" type="vector" simplifyDrawingTol="1" hasScaleBasedVisibilityFlag="0" simplifyLocal="1" maxScale="0" simplifyMaxScale="1">
      <extent>
        <xmin>58356.25057576912513468</xmin>
        <ymin>176062.54183909948915243</ymin>
        <xmax>66236.00066302371851634</xmax>
        <ymax>186876.7839431082829833</ymax>
      </extent>

And several other lines of code like it. Is it possible to make every project file automatically zoom to extents by changing something in the code? Or am I looking in the wrong place?

Edit: QGIS project is opened outside QGIS, on which QGIS starts. (I hope I said this correctly...)


texts: Which are the steps the union functon goes through?

texts: While skimming the union function I could see it works with intersecting bounding boxes, layer features and attributes. But how, exactly does this function work; which are the lines that check for intersection, compose new features, adds the attributes to those features, etc?


texts: Loading geopackage into mapView in Java?

texts:
Is there any way to load a geopackage file into any mapView in Java?

Is there any library to do this?


texts:
I can't find how to customize syntax highlighting for mapfile (*.map files) in Notepad++. I actually have add this syntax file recommended in this post. BUT the syntax highlighting is not good, probably because built for an old version of MapServer (see image below where plenty of key-words are not colored properly). I therefore have 2 questions:

How can I customize the syntax colorization? I can't access to the 'Mapfile' language via Settings > Style Configurator > Default theme (not in the languages list).
How is it possible to make Notepad++ recognize *.map files automatically. When I open such a file, it is colored like a txt file and the extension *.map is not recognized.


texts: Change connection details in QGIS 3 (Polygon to MultiPolygon)

texts:
I've got a qgis project with a lot of layers. In the past they were established as postgis-polygon-layers/-tables.

I updated my postgres/postgis-database recently.

Now the updated database stores the tables as MultiPolygon tables but QGIS does not recognize the modified type with the old connections. The change is only recognized if the tables are imported new into the project. It's the same with using qlr's with "Polygon"-connection-details.

I'd like to use the MultiPolygon-format to have less trouble cutting polygons into pieces.... So it would be nice to change the data-source-connection from Polygon to MultiPolygon.

As there are more than a hundered layers with relations to other layers in the project it would be a huge effort to replace the "old" Polygon layers with "new" MultiPolygon-layers and establish all the relataions, symbology ... manually

Is it possible to change the connection-details in a efficient, secure way (like python script / editing the qgs-file (replace all "Polygon" text with "MultiPolygon")) ?

I tried to change every layer-child of the layertree with the following code, but the connection is not changed.

If you run the code a second time the code does not enter the if-clause  - so I thing some changes are applied but not enough.

What do I have to do to apply my changes in the layer-tree ?

from PyQt5.QtCore import *
from qgis.gui import QgsMessageBar

def editConnectionDeeperInstance(laTrGr):
    pos = 0
    for child in laTrGr.children():
        if isinstance(child, QgsLayerTreeGroup):
            editConnectionDeeperInstance(child)
        elif isinstance(child, QgsLayerTreeLayer) and child.layer():
            l=child.layer()
            provider = l.dataProvider()
            if provider.name() == 'postgres':
                uri = provider.dataSourceUri()
                for t in tablesToEdit:
                    if re.search(t,uri) and re.search("type=Polygon",uri):
                        print('changed:',child.name(),'-', t)
                        new_uri=x = re.sub("type=Polygon", "type=MultiPolygon", uri)
                        # print(new_uri)
                        provider.setDataSourceUri(new_uri)
                        provider.reloadData()

tablesToEdit=["""public"."table1""","""public"."table2""" ]
laTrRo = QgsProject.instance().layerTreeRoot()
editConnectionDeeperInstance(laTrRo)


texts: Use WMTS as a mask to clip a polygon in QGIS

texts:
I have loaded a WMTS layer. But when I carry out some raster geo-algorithms nothing shows up in the drop down table (i.e. I can't use the raster layer as an input object).

How do I make it possible to use it as a mask to clip polygon?


texts: Joining 2D polylines and DEM into 3D polylines using QGIS/GRASS/SAGA?

texts:
I have a road network in 2D/PostGIS and a DEM. I want to turn the road network into 3D by retrieving the Z-value from the DEM.

Both files are quiet large so an efficient process is needed. I'd prefer to use QGIS or any other open source tool.


texts: Debugging ERROR 000728: Field *fieldname* does not exist within table even though field exists?

texts:
I do some script for a shapefile's  table.

And this script works for 2 fields and do not for other ones, because of:

ERROR 000728: Field FID_158_01 does not exist within table

But this field is fine and EXIST!

Script executed in field calculator.
I take field names from Iterator (Iterate MultiValue).

Expression:

a(!%field_name%!)

Codeblock:

def a(f):
  if f==0:
    z=1
    return z
  elif f==(-1):
    z=0
    return z
  elif f==(1):
    z=1
    return z

Crashlog:

Executing (EXTENT redo math): Model13 H:\TEST\CHECK.shp FID_154_01;FID_153_01;FID_158_01;FID_155_01;FID_153_02
Start Time: Fri Sep 28 14:12:22 2012
Executing (Iterate Multivalue): IterateMultivalue FID_154_01;FID_153_01;FID_158_01;FID_155_01;FID_153_02
Start Time: Fri Sep 28 14:12:23 2012
Succeeded at Fri Sep 28 14:12:23 2012 (Elapsed Time: 0,00 seconds)
Executing (Calculate Field): CalculateField H:\TEST\CHECK.shp FID_154_01 a(!FID_154_01!) PYTHON_9.3 "def a(f):\n  if f==0:\n    z=1\n    return z\n  elif f==(-1):\n    z=0\n    return z\n  elif f==(1):\n    z=1\n    return z"
Start Time: Fri Sep 28 14:12:23 2012
Succeeded at Fri Sep 28 14:12:23 2012 (Elapsed Time: 0,00 seconds)
Executing (Iterate Multivalue): IterateMultivalue FID_154_01;FID_153_01;FID_158_01;FID_155_01;FID_153_02
Start Time: Fri Sep 28 14:12:24 2012
Succeeded at Fri Sep 28 14:12:24 2012 (Elapsed Time: 0,00 seconds)
Executing (Calculate Field): CalculateField H:\TEST\CHECK.shp FID_153_01 a(!FID_153_01!) PYTHON_9.3 "def a(f):\n  if f==0:\n    z=1\n    return z\n  elif f==(-1):\n    z=0\n    return z\n  elif f==(1):\n    z=1\n    return z"
Start Time: Fri Sep 28 14:12:24 2012
Succeeded at Fri Sep 28 14:12:24 2012 (Elapsed Time: 0,00 seconds)
Executing (Iterate Multivalue): IterateMultivalue FID_154_01;FID_153_01;FID_158_01;FID_155_01;FID_153_02
Start Time: Fri Sep 28 14:12:24 2012
Succeeded at Fri Sep 28 14:12:24 2012 (Elapsed Time: 0,00 seconds)
Executing (Calculate Field): CalculateField H:\TEST\CHECK.shp FID_158_01 a(!FID_158_01!) PYTHON_9.3 "def a(f):\n  if f==0:\n    z=1\n    return z\n  elif f==(-1):\n    z=0\n    return z\n  elif f==(1):\n    z=1\n    return z"
Start Time: Fri Sep 28 14:12:25 2012
Failed to execute. Parameters are not valid.
ERROR 000728: Field FID_158_01 does not exist within table
Failed to execute (Calculate Field).
Failed at Fri Sep 28 14:12:25 2012 (Elapsed Time: 0,00 seconds)
Failed to execute (EXTENT redo math).
Failed at Fri Sep 28 14:12:25 2012 (Elapsed Time: 3,00 seconds)

Any ideas?

When i do all the steps in one model it gives me an error.

But when i just make union output file then stop the model and RUN THE SECOND part with field calculations -> error dissapeared and all works fine!!!


texts: Distinguish between filled polygon and boundary

texts:
I am trying to convert the raster images to vector and get the inner filled polygon. But the filled polygon and its boundary are wrongly detected for some cases. Below is the code where I tried to distinguish using area. How to distinguish them?
drv = ogr.GetDriverByName("ESRI Shapefile")
dataSource = drv.Open('test/raster_to_vector.shp', 1)

#print (dataSource)

layer = dataSource.GetLayer()

list1=[]

list2=[]

for feature in layer:
    geom = feature.GetGeometryRef()

    dn1 = feature.GetField('DN')

    area = geom.GetArea()

    if (area > 10 and area < 3000):

        list1.append(dn1)

        list2.append(area)

    else:

        layer.DeleteFeature(feature.GetFID())


texts: How to programmatically change the timeExtent of a Map & update the timeSlider?

texts:
I am developing with the ArcGIS JavaScript web API and ArcGIS server 10.0.

I have a time-aware layer in in my map, and also a timeslider. The full time extent of the map is a year, with data every month.

I want that when the map loads, the data is shown for that particular month, and the slider is also on that month.

I can programatically set the timeExtent on the map, using map.setTimeExtent(), but the timeslider does not auto-update to this date.

How should the timeExtent to be set, so that the timeSlider auto-updates along with it?

To see an reproducible example, go to the sample: http://help.arcgis.com/en/webapi/javascript/arcgis/samples/time_sliderwithdynamiclayer/index.html, and open up your firebug/developer console, and enter

map.setTimeExtent(map.timeSlider.fullTimeExtent);

The map will show you data for the entire period, but the time slider does not update.


texts:
Ok, so with ArcGIS 10, ESRI has depricated 'Address Standardizing' in the new locators, I know you can still use the 9.3.1 locators to do this; but I heavily suspect that after the 10.1 version we will see this available less and less.The FGDC Addressing Standard places a heavy value on good addresses having the Prefix, Prefix-Direction, Street-Name, Suffix, Suffic-Direction values for a good clean form. And my agency has a need for the same.Are there other tools out there that people are using to take structured addresses and parse them out fairly-reliably so they can be stored? Are any of them free?

*Update:*I came across this on CodePlex (US Address Parser) and it has given me a really good jumping off point. Its really nice in that its managed code that I can use all over the place and so far it has had a very high sucess rate at breaking down parts to be reviewed line by line and part by part.


texts:
I have imported several Shapefiles into PostgreSQL+PostGIS database. However, it seems that the database doesn't support metadata. I can't find a way to visit or edit metadata for the datasets. I don't know how to add attributes for the dataset (like author, time of creation, time of last edit, creator's affiliations, telephone number, etc.).

Is it possible to add or edit these data into PostGIS database?

EDIT: This question should remain open (i.e. no answer should be accepted).
EDIT2: I strongly oppose the opinion that metadata should be handled by the database designer or maintainer. Metadata is not arbitrarily designed. It has standards, which won't change in a specified occasion. Either the GIS platform or the DBMS should help end users implement metadata handling mechanics, including database design following given national or industry standards, reading and writing metadata in a user-friendly interface, and extensibility to customize fields. Even though people are still allowed to design a metadata table all from scratch, they simply tend to use predefined metadata schemas (which will be compatible to national and industry standards, with less or no errors than designing a table from scratch).


texts:
Nothing happens when I click the Georeferencer menu in QGIS 2.18.13. No window appears (menu Raster - Georeferencer - Georeferencer). I last used the georeferencer a couple of weeks ago, most likely before the last bug fix update.

I tried fixing the problem by renaming my .qgis2 directory, essentially resetting QGIS to "factory settings". This did not fix the problem. Nothing informative appears in log. It seems the georeferencer doesn't even to try to start.

My set-up: QGIS 2.18.13 (QGIS code revision 8cba0bb9eb) on Windows 10 64-bit, installed using the OSGeo4W installer.

Anyone else having issues with the Georeferencer in the latest QGIS 2.18 release?


texts: Automatically creating map for each column of attribute table in QGIS?

texts:
I have an atribute table with yearly data for 10 regions, is there a way how to automatically create a set of maps (one for each year) with fixed scale?

I tried to search the web but perhaps I didnt ask the correct question

Picture - how I would like the result to look like


texts: ArcGIS 10: How to clip point buffer to features in another layer?

texts:
I'm on a steep learnig curve here, so aplogies if my question is a one-banana problem.

In ArcGIS 10 I have a contour layer that was generated from an ESR ASCII raster file (batyhmetry of the South Pacific). I also have a point layer (indicating certain seabird breeding locations). For each point I produced 20km buffer zones. So far so good.

What I need to do now, is clip the circular buffer zones so that only the buffer regions between two contour lines remain (in my the case 0m and 150m isolines).

Is there a way in ArcGIS to clip/remove the buffer zones that are outside the 0-150m contours?


texts: How can I read in Qt listWidget functionality as Qgis the raster calculator?

texts:
I want to create a module in Qgis, my raster calculator, I'm stuck in the moment as it is possible to read raster files, just like in Qgis.
If we open the raster calculator Qgis, we see in the raster channel our raster files with which we can perform mathematical operations

In my code it reads the raster channel as plain text

def run(self):
    """Run method that performs all the real work"""
    layer_list = []
    layers = self.iface.legendInterface().layers()
    for layer in layers:
         layer_list.append(layer.name())
         self.dlg.listWidget.addItems(layer_list)

    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        pass

How now to specify what he would perceive in the listWidget, how exactly the raster (this is necessary so that later I could add the functions Gdal and read the channels of the raster)


texts: How to Change the Default Displaying Order of Esri Basemaps in basemapGallery?

texts:
I just started using ArcGIS JS APIs. In my app, I want to combine a few ESRI basemaps with my own basemap service. I am using basemapGallery dijit to change basemaps. By default, basemapGallery shows the imagery map. I want to show OpenStreetMap as the default. How can I make this happen?

The code in my app is simple and directly borrowed from ArcGIS samples. So I did not put it here.

Thanks!


texts: Calculating stats for single band raster using PyQGIS?

texts:
In the QGIS script editor, how do I calculate the main stats (minimum, maximum, mean, median, standard deviation,....) for a single band raster?

I am trying to make a tool for extracting a MASK (0-1) to cover clouds and cloud shadows for Landsat imagery. The algorithm is based on the 'MEAN' and 'Standard Deviation' of Bands BLUE and NIR. To do this I wrote the following script using the 'script editor' within 'Script' tool into the Processing Toolbox of the QGIS. Indeed the script works well, but I would like to automatically calculate 'MEAN' and 'Standard Deviation' instead of entering the values manually.

 ##Shadow-Cloudless MASK=name
 ##BLUE=raster
 ##NIR=raster
 ##MeanBLUE=number 0.0
 ##StDevBLUE=number 0.0
 ##MeanNIR=number 0.0
 ##StDevNIR=number 0.0
 ##StDevDivisionFactorBLUE=number 1.0
 ##StDevDivisionFactorNIR=number 1.0
 ##MASK=output raster

 from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
 import numpy as np

 #Get layer object
 layer1 = processing.getObject(BLUE)

 #Get layer object
 layer2 = processing.getObject(NIR)

 #Get number
 number1 = MeanBLUE

 #Get number
 number2 = StDevBLUE

 #Get number
 number3 = MeanNIR

 #Get number
 number4 = StDevNIR

 #Get number
 number5 = StDevDivisionFactorBLUE

 #Get number
 number6 = StDevDivisionFactorNIR

 def mask (BLUE,NIR,MeanBLUE,StDevBLUE,MeanNIR,StDevNIR,StDevDivisionFactorBLUE,StDevDivisionFactorNIR,output):
    entries=[]

    #define raster 1 ("BLUE")
    raster1=QgsRasterCalculatorEntry()
    raster1.ref='BLUE@1'
    raster1.raster=BLUE
    raster1.bandNumber=1
    entries.append(raster1)

    #define raster 2 ("NIR")
    raster2=QgsRasterCalculatorEntry()
    raster2.ref='NIR@1'
    raster2.raster=NIR
    raster2.bandNumber=1
    entries.append(raster2)

    number1='MeanBLUE'
    number2='StDevBLUE'
    number3='StDevDivisionFactorBLUE'
    number4='MeanNIR'
    number5='StDevNIR'
    number6='StDevDivisionFactorNIR'

    #MASK Processing
    calc=QgsRasterCalculator('((("BLUE@1"<('+str(MeanBLUE)+'+('+str(StDevBLUE)+'/'+str(StDevDivisionFactorBLUE)+')))+("NIR@1">('+str(MeanNIR)+'-('+str(StDevNIR)+'/'+str(StDevDivisionFactorNIR)+'))))=' + str(2) + ')',output,'GTiff',NIR.extent(),NIR.width(),NIR.height(),entries)
    calc.processCalculation()

mask(layer1,layer2,number1,number2,number3,number4,number5,number6,MASK)


texts: GeoPandas transform latitude/longitude coordinates to UTM coordinates

texts:
I have a shapefile which I read with GeoPandas, this shapefile has lat/long coordinates and I want to transform it to UTM coordinates.
I tried this:

shapefile = gpd.read_file('example.shp') #CRS = epsg:4326
shapefile.to_crs("+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs",inplace=True)

But it doesn't transform my coordinates correctly. I think it is because I have to explicitly put in the code the UTM zone with the number and letter which can vary depending on the shapefile I read.

So for calculating the letter and number I tried this:

import utm
utmZones = set()
for p in shapefile.geometry:
    utmZones.add(utm.from_latlon(p.centroid.y,p.centroid.x)[2:4])
print(utmZones)

    {(20, 'H'), (20, 'J'), (21, 'H')}

If the polygons fall in three different adjacent UTM zones is there any problem if I choose the UTM zone which is in the middle or should I create three different shapefiles one for each UTM zone?


texts:
I am thinking about a new QGIS plugin I'm going to do in near future.

So, I have some things, that could be hard for me. One of them is to making the plugin icon active when I click it.

For example, when I click on the "identify" tool in QGIS, it becomes active, but when I press the "pan" tool, the pan tool becomes active and the identify tool is no longer active.

I want do the same to my plugin icon.

Can you explain to me a little bit about how I'm supposed to do that?


texts: Analysing movement using cost surfaces with ArcGIS Path Distance tool?

texts:
I would like to analyse hypothetical movement (on foot) in a landscape as based on energy expenditure, but I have run into some trouble that I am hoping that you could help me with. I have tried to do this using ArcGIS’s Path Distance-tool in Spatial Analyst using cost surfaces I have created, but they result is not what I would have expected.

This is what my elevation surface looks like (downloaded from ASTER GDEM):

Based on the elevation data I created a cost surface that is supposed to contain energy expenditure (metabolic rate in Watts) per map unit (m). For this I used this formula:
M = 1.5W + 2.0 (W + L) (L / W)2 + N (W + L) (1.5V2 + 0.35V * abs(G + 6))

Or put in Raster Calculator terms:
(1.5 * 60) + (2.0 * (60 + 3) * Square((3 / 60))) + (1.2 * (60 + 3) * (Square((1.5 * "movementspeed")) + (0.35 * "movementspeed") * Abs(("slopeinpercent" + 6))))

Where M is the metabolic rate in Watts, W is the weight of the modelled individual, L is the individual’s carried load, N is a factor that describes ease of movement in the terrain (for testing purposes set to 1.2), V is the individual’s movement speed and G is the slope in percent. This created a surface with values ranging between 90 and 25000, with the majority of values between 90 and 1000 (which seems about right, the absurdly high values are most likely a result of flawed slope values, which could easily be fixed).

Movement speed was calculated using this formula:
V = 6e^(-3.5 * |s + 0.05| where s is the slope in degrees.

Or put in Raster Calculator terms:
6 * Exp( - 3.5 * Abs(Tan("slopeindegrees") + 0.05))
This created a surface with values between 0 and 5.9 km/h, which seems about right and is consistent with what I expected.

Now, these surfaces were used as input in the Path Distance-tool; the DEM as the input surface raster (i.e. in_surface_raster), the surface with energy expenditure as a cost raster and the DEM as the vertical raster to allow the tool to calculate whether or not the modelled individual is moving up or down a slope. For testing purposes, two points in the north-western and south-eastern corners of the DEM were used as source data (i.e. in_source_data). The output was this (red are unintuitively the lowest values and blue the highest):

My interpretation of the output is that it pretty much ignores differences in elevation, and the differences in value are simply related to differences in distance. I would have expected the surface to follow the flatter areas in the western part of the region and avoided the mountainous eastern parts, which it clearly doesn’t do. But, I’m still rather new to these types of analyses, and would appreciate the interpretations of others.
So, is anyone able to point out any flaws in my methodology/formulas that could cause the strange output? Or, is the output expected and I’m simply misunderstanding what I should expect from a path distance analysis?


texts: Multiple response from nominatim service call (Geocoding) in openstreetMap

texts:
I am using OpenStreeMap with Google as layer, and for geocoding , nominatim service.
Below, is my code.

$.ajax({
url: "http://nominatim.openstreetmap.org/search?q=" + address + &format=json&polygon=0&addressdetails=1",
type: 'POST',
//data: { 'htmlViewType': htmlViewType },
complete: function () { },
success: function (result) {
    var obj = JSON.parse(result);
    if (obj.length > 0) {
        document.getElementById('LblError').innerHTML = obj[0].display_name;
        var markerslonLat = new OpenLayers.LonLat(obj[0].lon, obj[0].lat).transform(WGS84, map.getProjectionObject(), 0);
        map.setCenter(markerslonLat, 10);
        //var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png', size, offset);
        var icon = new OpenLayers.Icon(' http://maps.google.com/mapfiles/ms/micons/blue.png', size, offset);
        markers.addMarker(new OpenLayers.Marker(markerslonLat, icon));
        map.addLayer(markers);
    }
    else {
        document.getElementById('LblError').innerHTML = 'no such address.';
    }
},
error: function (error) {
    document.getElementById('LblError').innerHTML = error;
}
});

Now, when i pass address = "Bharuch" in above URL , it returns 2 result in array, and when address =  "Surat", it returns 5.

so, my question is : which result value should i take from that array result as original Result ?

Right , now i am taking result from 0 index, but in few case it don't show my actual address.
Thanks.


texts:
I was looking for some publicly available ground data classification datasets, specifically, locations in e.g. latitude/longitude with a description of the land cover/land use, and a reference e.g. Landsat image. The aim of this is to have a sample classification dataset with which to try out different classification algorithms.

And yes, I could do this by hand, selecting random pixels from a landcover dataset, but I'd like to at least have some sort of ground truth to compare with and I'd rather use some freely available data that people have used for other things.


texts: How to create a WCS GetCapabilitiesType using Geotools?

texts:
I want to get the Capabilities of a WCS using Geotools (16.2) net.opengis.wcs20.GetCapabilitiesType. I found a WCS 2.0 test class on GitHub but was not able to use the code in the same way because this part:

GetCapabilitiesType caps = (GetCapabilitiesType) parser.parse(getClass()
                .getResourceAsStream(capRequestPath));

throws

  Exception in thread "main" java.net.MalformedURLException  [...] Caused by: java.lang.NullPointerException

I wrote this method to overcome the exception:

GetCapabilitiesType getCapabilitiesType(String capRequestPath)
        throws IOException, SAXException, ParserConfigurationException {
    Parser parser = new Parser(new WCSConfiguration());
    FileInputStream fis = new FileInputStream(capRequestPath);
    GetCapabilitiesType capabilitiesType = (GetCapabilitiesType) parser.parse(fis);
    return capabilitiesType;
}

This throws

  Exception in thread "main" java.lang.ClassCastException:
  java.util.HashMap cannot be cast to
  net.opengis.wcs20.GetCapabilitiesType

which I wasn't able to resolve. I also found a WCS 1.1 test class on GitHub.

This approach works:

HashMap<String, Object> getCapabilitiesMap(String capRequestPath)
        throws IOException, SAXException, ParserConfigurationException {
    Parser parser = new Parser(new WCSConfiguration());
    FileInputStream fis = new FileInputStream(capRequestPath);
    HashMap<String, Object> capabilitiesMap = (HashMap<String, Object>) parser.parse(fis);
    return capabilitiesMap;
}

but this way I do not have a net.opengis.wcs20.GetCapabilitiesType object.

Then I found the class net.opengis.wcs20.impl.Wcs20FactoryImpl on GitHub which has a createCapabilitiesType() method to create a net.opengis.wcs20.GetCapabilitiesType object. I created my GetCapabilitiesType like this:

String wcsURL = "https://www.wcs.nrw.de/geobasis/wcs_nw_dgm";
Wcs20Factory factory = Wcs20FactoryImpl.init();
GetCapabilitiesType wcsGetCapabilities = factory.createGetCapabilitiesType();
wcsGetCapabilities.setBaseUrl(wcsURL);

But when I want to use getters like:

List<String> versions = wcsGetCapabilities.getAcceptVersions().getVersion();

or

List<String> sections = wcsGetCapabilities.getSections().getSection();

they throw

  java.lang.NullPointerException

What is the best way to get a net.opengis.wcs20.GetCapabilitiesType object?


texts:
I want to iterate through a layer that has about 500 records. In my code, I'm using a definition query inside a for loop to iterate through a test list:

listx = [1,2,3,4,5]

My question (probably an elementary one), is how to just call out all records instead of manually typing all 500 object ID's? I want to iterate through all Object IDs, while ingesting each ID into the definition query. I thought trying [:] or [*] or [0:500] would work, but they don't.

import arcpy
from arcpy import mp

p  = arcpy.mp.ArcGISProject("CURRENT")
m = p.listMaps("Map")[0]

lyr1 = m.listLayers("No_Other")[0]
lyr2 = m.listLayers("No_Admin")[0]

#obid = "OBJECTID in (409, 764, 114, 776)"
obid = "OBJECTID = {}"
listx = [114, 409, 764, 776]

for oid in listx:
    lyr2.definitionQuery = obid.format(oid)
    arcpy.SelectLayerByLocation_management(lyr1, 'WITHIN_A_DISTANCE', lyr2, "150 Miles", 'ADD_TO_SELECTION', 'NOT_INVERT')
    lyr2.definitionQuert = None
    print(lyr2.definitionQuery)


texts:
I am working with Leaflet and EasyButton and have several EasyButtons on the map. I would like to move them outside the map, but have difficulty to achieve this. Seems EasyButton does have any function for capturing the name of the buttons that it holds.

Something similar to this, but with easyButton.
Placing controls outside map container with Leaflet?

var ButtonPower_Line = L.easyButton({
    id: 'ButtonPowerLine',
    states: [{
        stateName: 'Show',
        icon: ' Off',
        title: 'Show',
        onClick: function (btn) {
           //do sth
        }
    }, {
        stateName: 'Hide',
        icon: ' On',
        title: 'Hide',
        onClick: function (btn) {
          //do sth
        }
    }]
});
ButtonPower_Line.addTo(map);


texts: Calculating area of intersecting polygons within severals overlapping buffers

texts:
I am working with..

an initial point layer (in blue) called : centroid_layer.shp
an initial polygon layer (in green) of woods called : wood_zone.shp

..and I need to get the area of intersected wood_zone within each buffer (in red, created from centroid_layer), knowing that almost all my buffers are overlapping each others. Thus I need to include the wood_zone polygones as many times as the number of intersecting buffers.

I believed it is not possible to do only with QGIS processing and/or modeler (???) so I would like to solve it through PostGIS/PostgreSQL (9.6 version).

I have seen in a similar posted question from Anna and I was trying to proceed with SS_Rebelious solution using his SQL script not in Spatialite but in PostGIS/PostgreSQL (with the Query tool of pqAdmin 4) as following :

CREATE TABLE buff AS
SELECT gid, ST_Buffer(geom, 1000,'quad_segs=100') FROM centroid_layer ;

and then :

UPDATE centroid_layer
SET areacolumn = (SELECT ST_Area(
                          (SELECT ST_Intersection(
                            (SELECT ST_Union(geom) FROM wood_zone),
                            (SELECT geom FROM buff WHERE centroid_layer.gid = buff.gid)
                                                 )

                           )
                                )
                  );

but I have a problem after the UPDATE...SET... step. Indeed my resulting areacolumn is completly filled with 0.0 values (see below).

So I was wondering if I am doing something wrong or if there would be a different SQL script to use with Postgres-PgAdmin ?

EDIT : I have tried the same SQL script with a sample group of three buffers (associated with 3 related points) and I was able to collect the area values within my column named "areacolumn". At this point the only difference I see is that I used a different primary key, I mean I did not used the "gid" (automatically created by postgres/postGIS) primary key.

But I still do not understand... is the use of the "gid" as primary key a problem with such a SQL script ?


texts: I'm trying to import a GDB folder, downloaded from https://data.geo.admin.ch/ch.swisstopo.geologie-felsoberflaeche_hoehenmodell/ into QGIS 3.8 but I get an errormessage that data source is invalid. I have tested with QGIS on MacOS and also on OpenSuse. I have also tried to use ogr2ogr to get shape files. But also here I failed. Unfortunately, I cannot find out whether this file is corrupt or what are the problems to impede the working with these data.


texts: Is there a maximun number of edges/vertices PGRouting can cope with?

texts:
We have imported OSM networks to route on. So our IDs sometimes grow abnormally big ;)

It seems that PGRouting is failing with IDs like 633719916...

We got this error message from postgresql:
"terminate called after throwing an instance of 'std::bad_alloc'".

Is there another way rather than renumbering our IDs?

Is there a "fix" planned from PGRouting?

Thanks for your help!


texts: Every time I export a new print composer (either image, svg, or PDF) in QGIS 2.18.1, I lose my composer

texts:
As the title suggests, whenever I attempt to make a new print composer (in either windows or Linux, I have tried both) with my map, it goes blank and the file is not exported.

When I attempt to export as a PDF, I get the message "Advanced composition effects such as blend modes or vector layer transparency are enabled in this project, which cannot be printed as vectors. Printing as a raster is recommended." Would this prevent my ability to make a PDF from the print composer?

If so, how can I rectify this?

Since I have a low reputation, if you need pictures I can provide them in the comments.


texts:
I am new to DBA and PostGIS. I am developing an application using PostgreSQL and PostGIS where nearby locations need to be queried.

I created a sample table as below:

CREATE TABLE foo (geog geography);

Inserting like:

INSERT INTO foo (geog)
  VALUES (ST_MakePoint(12.234,23.22));

But in some tutorials,insertion is like:

INSERT INTO foo (geog) VALUES (ST_SetSRID(ST_MakePoint(-122.079513,45.607703), 4326));

I can't clearly understand what ST_SetSRID does.
What does it do?
Do I need to use it If I am going to find nearby places, etc?
If I should use it, what value should I go with in the place of 4326 ?


texts: Creating 3D polyline just below another polyline using ArcGIS Desktop?

texts:
I have a geological setting where I would like to add a lithology to an existing map. I already have the 3D polylines from my geological map and the DEM I now want to add another polyline between two outcrops that occurs like 20 m below one polyline and follows the topography.

Picture one shows the 3D model and two the 2D one. The black line is what i would like to create automated with ArcGIS Desktop.

Anyone an idea how to do that?


texts: How to create a dense pointcloud of trees in ArcGIS?

texts:
I'm having a large dataset of forest polygons and would like to scatter individual tree locations over these areas. Here is how I have done this so far:

In "Create Random Points" number of points has been set to 999999 and the minimum allowed distance (linear) to 5 meters.

Why am I running the "create random points" 4 times in that model?
The problem is that when running it once in ArcGIS, the points created are way too sparse. Running it four times and then merging all positions that are closer together than 1.4 meters (which essentially comes out to ~5 meters as tests have shown) gives me a much denser point cloud.

Any ideas how to change this approach into something faster with improved performance?


texts: Navigating real time with USB GPS usiang QGIS

texts:
I'm trying to navigate real time with an USB GPS over a shapefile map, but I can't find the way to do this. I'm using QGIS 1.9.0 desktop version with a Garmin Etrex Legend HCx GPS.

In fisrt place, QGIS doesn't recognise my GPS port. In addition, I cannot find in QGIS the option to navigate real time.

Do you have any idea to help me with this?


texts: Replace String With Number In List, Google Earth Engine?

texts:
How can I replace '-infinity' with 0 number using a function in Google Earth Engine? I'm gonna to replace string with 0 value in the list automatically.
https://code.earthengine.google.com/1bd2127a2d5f0994dfd593be190edd60
var x  = ee.List([0,1,2,3,4,5,6,'-infinity','-infinity']);

var y = ee.List([1,2,3,4,5,6,7,8]);

print('xvalues',x);
print('yvalues',y);


texts: How to calculate the shortest path of adjacent points, from an origin to destination points? (ArcGIS Pro/ArcMap)

texts:
For a project I am trying calculate the shortest path (marked by points, not lines) connecting a point of origin to a destination point (different feature layer) and automate this for several segments/paths.
The output should be the respective points as a new feature layer/shapefile and not lines, just the points with sequential number (from origin to destination).
I am familiar and have access to ArcGIS Pro/ArcMap, ArcPy or ModelBuilder/Geoprocessing tools.
Below I more detailed example:

I already considered using the near function from the ArcGIS geoprocessing options to calculate the distances between the Origin and closest points. However, simply following the closest distance of respective (blue) points, does not allow to construct the path that I marked in red. The yellow line is just for illustration, and should not have a function in the analysis
I would have to iterate/loop a near function for every consecutive blue point, while also simultaneously keep the angle (in line with the destination point) into account.
For example, if I would draw a line from nearest to consecutive nearest point, after the first blue point, the line would go down, instead of in line with the destination point. Which ultimately would be the shortest path.


texts: Legend Limiter will not install to my ArcView 10 - any suggestions?

texts:
In an older version of ArcView, I had an add-on called Legend Limiter for my software.  It would limit the legend categories to ONLY those shown in your map extent.  Say for example, I have a Geology shapefile containing over 100 categories of Geology, but the extent I'm showing on my map only includes 15 categories, this tool will limit my legend to reflect ONLY the categories in my map extent.  It is a great little time saver.

I installed the version found HERE.  When I try to run the .exe file however, I get the prompt shown below:

I tried googling the error message number, and the information does not resolve my issue.  I did find a fix for the newer ArcView HERE.  I follow the instructions from the forum to the letter, and when I add the .dll file, it's not recognized (prompt below).

I am officially out of ideas.  I just can't get it to work.  Any suggestions would be greatly appreciated.  Thanks all.


texts: Is there an ADO.NET 2.0 driver for the ESRI file geodatabase?

texts:
I am trying to make a connection to a file geodatabase to issue some SQL queries.

I don't need any spatial functionality, we're just trying to get data out of tables.  The application making the queries can use ADO.NET 2.0 drivers, which implement the DbProviderFactory classes to connect to various database and data sources.

I have found this: https://filegdbdotnet.codeplex.com/ but it's ADO 1.1 and doesn't implement the factory classes.

There's an Ole DB provider (ESRI.GeoDB.OleDB.1) but I think that requires an ESRI license to activate.  I'd like to avoid using that if possible since we can't be sure that the target machine running the application will have that.

Is there an ADO 2.0 driver for the file geodatabase?


texts: Change from unkown to a known coordinate system for a self made polygon

texts:
I have made a map with self made polygons. At the beginning I was so stupid I forgot to add a coordinate system to the new shapefiles. My map is finished now, but I can't export it because the Prj file is missing now.

When I'm trying to change it at the shapefile properties (XY coordinates) the polygon disappears or it moves to someplace in Africa (it is a map from Europe). Does somebody know a solutions for this?


texts:
I am trying to make a custom function in the QGIS expression editor that takes three arguments. Below is the tested standalone function, which I would like to integrate:

import math

def qfull_circle(DN, kb, I):
    #Input: DN in mm, kb in mm, I in ppm
    #Output: Q in l/s
    viscosity = 1.31 * math.pow(10,-6)
    roughness = kb / 1000
    Acircle = (DN * DN * math.pi / 4) * math.pow(10,-6)
    Xroot = math.pow((2 * 9.81 * DN * math.pow(10,-3) * I * math.pow(10,-3)),0.5)
    fraction1 = (2.51 * viscosity) / (DN * math.pow(10,-3) * Xroot)
    fraction2 = roughness / (3.71 * DN * math.pow(10,-3))
    return Acircle * (-2 * math.log(fraction1 + fraction2, 10) * Xroot) * math.pow(10,3)

Unfortunately, I cannot make it work. Any assistance on how to define it within the function editor is highly welcome.


texts: Line String Self-Intersection points in PostGIS

texts:
I'm trying to get all the points from a Self-Intersection Line String but not find out the self intersection point in the below attached image.

The intersection point touches on the same feature how to find out the self intersection point in postGIS?


texts:
I am working to upload a TIFF file into GeoServer so that it can be viewed in layer preview. I exported the TIFF from DroneDeploy, it has an associated geographic location. I was not able to get the GeoTIFF Store to recognize the URL of my file...kept getting errors. I tried using the raster2pgsql extension to upload into pgAdmin and then use a PostGIS store. I was able to upload, however, I get an error in layer preview, "Error rendering coverage on the fast path java.lang.NullPointerException"

I am semi-new to the open-source GIS world.


texts: Exporting rain data to CSV file from Google Earth Engine

texts:
I would like to get the CSV file of rain daily data from GSMAP for a specified polygon (mean in geometry).
I succeed in getting into it on a chart, but now I want to download the data and I'm having some trouble
Could you help me?
//Add Dates range, with las date available
var iniDate = ee.Date('2018-03-01')
var endDate = ee.Date(ee.List(GSMAP.get('date_range')).get(1)) //last date

print('lastDay available', ee.Date(ee.List(GSMAP.get('date_range')).get(1)))

//Add ImageCollection
var filterGSMAP = GSMAP.filterDate(iniDate, endDate).select(['hourlyPrecipRateGC']);

//Initiate the date list with a day less than the availability,
//as the last day available does not contain images:
var difdate = endDate.advance(-1, 'day').difference(iniDate, 'day'); //calculate steps/days in between
var createList = ee.List.sequence(0, difdate) //make list with the needed number of entries
var listdates = createList.map(function(day){ //change list into list of dates
return iniDate.advance(day, 'day')})
print(listdates)

//Summarize hourly precipitation data to daily precipitation data
//Then filter the image collection using the filterDate, which is (I think) easier than filtering
//on year and doy consecutively:
var gsmapImageCollection  = ee.ImageCollection.fromImages(listdates.map(function(summarize_day){
var filterCol = filterGSMAP.filterDate(ee.Date(summarize_day), ee.Date(summarize_day).advance(1, 'day'))
return filterCol.sum().select(['hourlyPrecipRateGC'], ['dailyPrecipRateGC'])
.copyProperties(filterCol.first()).setMulti({
Date: ee.Date(summarize_day), 'system:time_start': ee.Date(summarize_day).millis()});}))
print(gsmapImageCollection)
print(ui.Chart.image.series(gsmapImageCollection, geometry, ee.Reducer.mean(), 5000)
.setChartType('ColumnChart'));

//Set a system:time_start property,
//so you can easily plot/filter etc the resulting image collection.
//Add the resulting image collection to the map
var VIS = {palette: ['1621a2','ffffff', '03ffff','13ff03','efff00','ffb103', 'ff2300'], min: 0, max: 50}
//Map.addLayer(gsmapImageCollection.first(), VIS, 'GSMap');

// Here is my doubt, how can i get the CSV file with two columns, day and acumulated rain.
//Export the time-series as a csv.
Export.table.toDrive({
  collection: gsmapImageCollection,
  description: 'rain',
  selectors: 'day, value',
  fileFormat: 'CSV'
});


texts: gdal2tiles.py multithreaded only works with base tileset

texts:
I'm using GDAL 2.4.0 on Ubuntu 18.04 Installed from the ubuntugis-unstable repository.

When I use the gdal2tiles.py with the '--processes=n' flag the multithreaded operation works fine for processing the base tileset. Then, generating the overview tiles is only single threaded. Is this the expected behavior, or is something amiss in my configuration?


texts: Avoiding fails from ArcObjects geoprocessing with .NET?

texts:
There are some nice features in ArcToolbox we can use, but for some reason, this is NOT working properly. It doesn't even throw me an error.

My software is running inside ArcMap, so no need to AoInitialize again, corret?

    public void Execute()
    {
        InitializeProduct();
        try
        {
            Geoprocessor gp = new Geoprocessor();
            gp.OverwriteOutput = true;

            FeatureToPoint featureToPoint = new FeatureToPoint();

            string outputPathName = CurrentWorkspace.PathName + "\\teste_centroide";

            featureToPoint.in_features = InputFeatureClass;
            featureToPoint.out_feature_class = outputPathName;
            featureToPoint.point_location = "INSIDE";

            IGeoProcessorResult result = (IGeoProcessorResult)gp.Execute(featureToPoint, null);

            if (result == null)
            {
                for (int i = 0; i <= gp.MessageCount - 1; i++)
                {
                    Console.WriteLine(gp.GetMessage(i));
                }
            }

            IGPUtilities gpUtils = new GPUtilitiesClass();
            this.OutputFeatureClass = gpUtils.OpenFeatureClassFromString(outputPathName);
        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.Message + "\r\n");
        }

This is a code example I'm having here. I generated the DataManagement tools assembly, but I could not find the file to sign it.

This code just gives me an error. is it because of the signing?

I've tried the other way too, using IVariantArray and calling from the tool name, without sucess. Is it just me or...?

Can anyone point me a "nicer" solution? I need to run several process that are already built in ArcToolbox that I really don't want to duplicate.


texts: Dropping OID column with ArcPy Table To Excel?

texts:
Using ArcPy in a script, I'm exporting a table to Excel (using Table To Table and Table to Excel).

Is there a way or workaround to get rid of the Object ID column in the export process?

I'd like to avoid using the Python Excel modules.


texts: How to create tile images for tiled maps using ArcGIS Online?

texts:
I am a newbie to ArcGIS Online. I have to create a temperature tile map based on my own custom data. I have found that a tile layer can be added using the following link:

http://doc.arcgis.com/en/arcgis-online/reference/tile-layers.htm

Adding a new ArcGIS tile layer uses a service request like

(http://{subDomain}.tile.opencyclemap.org/cycle/{level}/{col}/{row}.png)

to get specific images from service for specific zoom levels, the question is how and where these images will be created based on our own data?

Is there any tool that creates such images for storage purpose and later to utilize using service call?

Can we create these images dynamically at run time based on our own data and how?
We want to create the end product like:
http://maps.google.com/gallery/details?id=zttWLnOPAlrs.kQBKSYw5ok5U&hl=en


texts: Loading and georeferencing historic maps using QGIS?

texts:
I have recently started working with QGIS.

I am unsure of how to go about inserting historical maps which have no referencing data.

I think it might be possible to insert them, then move and scale them using a modern mapping layer for reference.


texts: Issue with gdalUtils vs. raster in gridding data in R

texts:
I have irregularly spaced data (XYZ) and I'm trying to grid it in R using 2 different libraries, gdalUtils and raster.  I'd prefer to use gdalUtils as I have more options, but the output looks off.  To use gdalUtils I create a CSV file with columns named x, y, and z.  I then create a VRT file so that the GDAL function gdal_grid can be used.  My code,

## load libraries
library(ncdf4); library(raster); library(gdalUtils)

## download file and open
file = 'http://ncss.hycom.org/thredds/ncss/grid/GLBa0.08/expt_90.9/2012?var=mld&north=47&west=277&east=320&south=22&disableProjSubset=on&horizStride=1&time=2012-01-01T00%3A00%3A00Z&vertCoord=&addLatLon=true&accept=netcdf'
saveFile = tempfile(, fileext = 'nc')
download.file(file, saveFile, mode = 'wb')
ncOpen = nc_open(saveFile)

## get data
x = ncvar_get(nc = ncOpen, varid = 'Longitude') - 360
y = ncvar_get(nc = ncOpen, varid = 'Latitude')
z = ncvar_get(nc = ncOpen, varid = 'mld')

## replace NA with a fillValue
z[is.na(z)] = -9999

## create a data.frame
df = data.frame(x = c(x), y = rev(c(y)), z = c(z))

## close .nc file
nc_close(ncOpen)

## create a temp directory
temporary_dir = getwd()

## base dir path
tempfname_base = file.path(temporary_dir, 'hycom')

## create a properly formatted CSV
tempfname_csv = paste0(tempfname_base, '.csv')
write.csv(df, file = tempfname_csv, row.names = FALSE)

## tiff file for saving
tempfname_tif = paste0(tempfname_base, '.tiff')

## create VRT file
tempfname_vrt = paste0(tempfname_base, '.vrt')
vrt_header = c(
    '<OGRVRTDataSource>',
    '\t<OGRVRTLayer name="hycom">',
    '\t<SrcDataSource>hycom.csv</SrcDataSource>',
    '\t<GeometryType>wkbPoint</GeometryType>',
    '\t<GeometryField encoding="PointFromColumns" x="x" y="y" z="z"/>',
    '\t</OGRVRTLayer>',
    '\t</OGRVRTDataSource>'
)
vrt_filecon = file(tempfname_vrt, 'w')
writeLines(vrt_header, con = vrt_filecon)
close(vrt_filecon)

## now run gdal_grid
ras = gdal_grid(src_datasource = tempfname_vrt, dst_filename = tempfname_tif, a = paste0('nearest:nodata=', -9999, ':radius1=0.08333333:radius2=0.08333333'), txe = c(-83, -40), tye = c(22, 47), outsize = c(538, 385), a_srs = '+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0', of = 'GTiff', ot = 'Float64', l = 'hycom', output_Raster = TRUE, verbose = FALSE)

Everything runs fine, but when you plot the output raster and overlay the U.S. coastline you can see there is an issue.

However, if I use rasterize from the raster library everything looks good.  Here's my code,

x = ncvar_get(nc = ncOpen, varid = 'Longitude') - 360
y = ncvar_get(nc = ncOpen, varid = 'Latitude')
z = ncvar_get(nc = ncOpen, varid = 'mld')
df = data.frame(x = c(x), y = c(y), z = c(z))
eRas = raster(ncol = 538, nrow = 385)
extent(eRas) = c(-83, -40, 22, 47)
proj4string(eRas) = '+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0'
ras = rasterize(df[,c('x', 'y')], eRas, df$z)

I can't figure out why gdal_grid is giving me an image that looks offset.  The input data is exactly the same in both cases.


texts:
Is there a straight method to export shapefile geometry (polygon) to WKT using pyshp module for python?

I was looking for something similar to OGR's .ExportToWkt() method, but it seems like there's nothing similar implemented in pyshp


texts: I am looking for an example of how to use the QgsRuleBasedRendererV2 to filter features of a layer based on an attribute value (0 or 1). Has anyone done this in Python or C++?


texts:
I am getting the below error during execution of archydro.
"system runtime interoperservices com exception (0x80004005): error HRESULT E_FAIL"

I tried all options available on internet to solve this but no luck.

What is the cause of these type of issues? Is this related to application or os ?
What is the perfect fix if it is caused by arcgis/hydro or os ?

I am using windows7, 64-bit, Arcgis10.


texts: QGIS project wise grant to Users created in PostgreSQL

texts:
I have created multiple projects IN QGIS  and saved them in PostgreSQL. It has created a table in Public schema as 'QGIS_Projects' and creates rows to save different projects from QGIS . I have given 'Select' grant to users for this project table. I am now unable to give grants to users based on the different QGIS projects :
Example: QGIS _PROJECT 1 -- Select grant only to user A
Example: QGIS _PROJECT 2 -- Select grant only to user B
How can I do this?


texts: How to add vertex to an existing polyline programatically?

texts:
I have a Polyline and would like to add new vertex to before start point or end point or at any part of the polyline using python in qgis.

How can i do that?


texts: I have downloaded the heron mc and placed in my local C Drive and tried running some of the example maps. But all of the templates that use getfeatureinfo (getfeatureinfopanel, etc) return no features found. By doing some research I think the issue is I need to setup a proxy in Apache and in the proxy.cgi file. However I am not clear on how to do this, nothing I am trying is working or allowing the maps to bring up the features.


texts:
I have the same problem as this How to set up QGIS print composer options from the python console?

I have a composer template called "Template", I can open it with the python console but not edit the map content.

I'm very new to python so its likley to be something simple Ive overlooked. The ultimate aim is to loop through locations (from a shapefile) and create a map for each one with unique labels. I realise people have done similar things before, but none of my attempts have got them to work so I'm starting from scratch.

I'm using QGIS 2.8 on Windows 7.

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtXml import *
import tempfile
import os

#test
i=1

cLayer = iface.mapCanvas().currentLayer()
expr = QgsExpression( "\"ID\"=i" )
it = cLayer.getFeatures( QgsFeatureRequest( expr ) )
ids = [i.id() for i in it]
cLayer.setSelectedFeatures( ids )

canvas = qgis.utils.iface.mapCanvas()
canvas.zoomToSelected()

myFile = os.path.join('Z:\Output', 'Template.qpt')
myTemplateFile = file(myFile, 'rt')
myTemplateContent = myTemplateFile.read()
myTemplateFile.close()
myDocument = QDomDocument()
myDocument.setContent(myTemplateContent, False)
(True, u'', 0, 0)
newcomp = iface.createNewComposer()
newcomp.composition().loadFromTemplate(myDocument)
True
composition = QgsComposition(canvas.mapSettings())
map_item = composition.getComposerItemById('map')
map_item.setMapCanvas(canvas)
**Traceback (most recent call last):
File "<input>", line 1, in <module>
AttributeError: 'NoneType' object has no attribute 'setMapCanvas'
map_item.zoomToExtent(canvas.extent())
Traceback (most recent call last):
File "<input>", line 1, in <module>
AttributeError: 'NoneType' object has no attribute 'zoomToExtent'***


texts: Are there any shapefiles available online of ZCTAs (Zip code tabulation area) by state? The US Census bureau has an excellent list of ZCTAs but they do not indicate which state each ZCTA is in. They have a set of files that describe relationships between ZCTA and other districts, but not states.


texts: Difference between planar and Euclidean distance?

texts:
I wanted a straight-line distance. When I was calculating near distance on Arc GIS, I was only given the option between planar and geodesic distance.

But people in my field keep referring to my straight line distance as Euclidean.

I understand that planar is the default for my coordinate system (I think), and that it ignores surface curvature (which I wanted to ignore).

I am just worried that I calculated the wrong type of distance, but Euclidean was never an option in Arc GIS.


texts: Any ideas to display SLD TextSymbolizer property in OpenLayers

texts:
I'm trying to display the SLD TextSymbolizer property for a vector layer in OpenLayers but I found that it's not supported.

Has anyone ever dealt with this issue?

I found these links but I have not ideas about how to implement them

http://trac.osgeo.org/openlayers/ticket/1542

http://trac.osgeo.org/openlayers/ticket/2816


texts: Open a NetCDF file with multiple time step on QGIS (using Python console)

texts:
How to load a NetCDF-file into QGIS using Python? explains well how to load a NetCDF with PyQgis. However, it did not help me to choose among all the available time steps.

Is there a way to import a specific time step (or all of them) when I am loading layers from a NetCDF file?


texts: GeoExplorer cannot upload shapefiles into PostGIS

texts:
I have installed OpenGeo 4.1 on Ubuntu Linux 12.04 following the official instructions.

I then followed the tutorial Load your data with GeoExplorer.

I was able to zip and upload a few shapefiles. The layers show at the right geographical extents. However, what I didn't get is the following, as said in  the tutorial

  When deployed alongside GeoServer and PostGIS, GeoExplorer will
  convert spatial files into PostGIS tables, and then serve that data
  through GeoServer. This process is seamless; all you need to do is
  supply the data.

The shapefiles I uploaded didn't become PostGIS tables as advertised. Instead, when inspecting the uploaded layers in GeoServer, I see that the data source under "Connection Parameters" for the layers are still something like this (a shapefile):

file:uploads/tmp6091327178585086343/nr_adm0/nr_adm0.shp

I am able to access PostreSQL from pgAdmin, as well as from command line with

psql -U postgres -h localhost

I was wondering if someone has had a similar problem and know where the issue might be with the link from GeoExplorer to PostGIS?


texts: Getting string attribute from a shapefile with GDAL Feature.GetFieldAsString(name) returns incomplete value

texts:
I read a shapefile (zipped) with GDAL library in C#. Shapefile contains .cpg file with "UTF-8" encoding specified.

When I call Feature.GetFieldAsString(attributeName) some of the attribute values get trimmed at the end, and some have broken symbols.

Converting shapefile to GeoJSON with ogr2ogr makes a GeoJSON file with all attributes as they are supposed to be.

For example, value returned by GetFieldAsString is "\"Тер�", when the real value, also seen in geojson is "\"Терем+\"".

This doesn't seem to be an issue with encoding, for the beginning of the string is read correctly. But I have no clue where to look for the solution for this issue.


texts:
I've done a lot spatial analysis in ArcGIS. Then I decided to start to learn QGIS. I started with DEM analysis. But I found differences in statistics of DEM. ArcGIS correctly recognizes the lowest and maximum values (331 m, 2557 m), but QGIS does not. Based on statistics in QGIS the lowest value is 17 m greater (348 m) and the largest value 126 m greater than values in ArcGIS.

I also see differenet values when adding DEM to the table of contents in QGIS. The lowest value is 398 m and the largest 2066 m. Can anybody help me how to calculate correct values?


texts:
I have a variogram exported from Supervisor which contains (in this case) the Nugget, Sill and Range for two Structures in four directions.

How can I use this as inputs into a plot?

Is there some way to calculate a "fake" covariance for certain distances so I can draw a line between XY points?


texts:
I am unable to display ArcGIS dynamic map in Android but same code works for Tile maps.

Error I am getting is

  Caused by: java.lang.UnsatisfiedLinkError: Couldn't load GLMapCore
  from loader dalvik.system.PathClassLoader[DexPathList[[zip file
  "/data/app/com.esri.arcgis.android.samples.helloworld-1.apk"],nativeLibraryDirectories=[/data/app-lib/com.esri.arcgis.android.samples.helloworld-1,
  /system/lib]]]: findLibrary returned null


texts: Jquery attr() funciton returning "undefined" for GetCapabilities Request

texts:
On my webpage where I am displaying my maps, I want to add a link to the metadata of the selected map.
I have specified the link in GeoServer in the MetaDataURL field, and should be able to retrieve it with a GET request and finding it in the XML(See below for XML and Jquery code)
But every time I test it I get "undefined"...which is very strange, it seems as if the attribute doesn't even exist. Does anyone have a clue why it returns "undefined"

$.ajax({
        type: "GET",
        url: "http://fast.openearth.eu/geoserver/wms?service=wms&version=1.1.1&request=GetCapabilities",
        dataType: "xml",
                       success: function(xml) {
                                               $(xml).find('Layer').each(function(){
                                                     if ($(this).children("Title").text() == layer.name) {
MetaDataURL = $(this).children("MetadataURL").children("OnlineResource").attr("xlink:href");
}});}});

<Layer queryable="1">
                                                                                <Name>Global:wcmc-010-mangroveusgs2011</Name>
<Title>Mangrove USGS2011</Title>
<Abstract>Giri C, Ochieng E, Tieszen LL, Zhu Z, Singh A, Loveland T, Masek J, Duke N (2011).&#13;
    Status and distribution of mangrove forests of the world using earth observation&#13;
    satellite data. Global Ecology and Biogeography 20: 154-159. URL:&#13;
    http://onlinelibrary.wiley.com/doi/10.1111/j.1466-8238.2010.00584.x/abstract ;&#13;
    http://data.unep-wcmc.org/datasets/4</Abstract>
<KeywordList>
<Keyword>wcmc-010-mangroveusgs2011</Keyword>
<BoundingBox CRS="CRS:84" minx="-180.537994384766" miny="-35.6819496154785" maxx="181.735992431641" maxy="30.3975315093994"/>
<BoundingBox CRS="EPSG:4326" minx="-35.6819496154785" miny="-180.537994384766" maxx="30.3975315093994" maxy="181.735992431641"/>
<MetadataURL type="ISO19115:2003">
<Format>text/plain</Format>
<OnlineResource xlink:type="simple" xlink:href="http://fast.openearth.eu/geonetwork/srv/eng/metadata.show?id=39&amp;currTab=simple"/>
</MetadataURL>


texts: EPSG:31254 to WGS84 Transformation Parameters

texts:
The official transformation parameters from EPSG:31254 to WGS84 are:

  towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 (EPSG:1618)

(defined here and used by ArcGIS MGI_To_WGS_1984_3 transformation method as well as in the EPSG:31254 definition used by QGis 2.18.1)

GeoServer 2.11.0 on the other hand uses the following transformation parameters in its EPSG:31254 definition:

  TOWGS84[601.705, 84.263, 485.227, 4.7354, -1.3145, -5.393, -2.3887].

If I provide GeoServer with an alternate version of EPSG:31254 by using the epsg_overrides.properties file - using the EPSG:1618 transformation parameters -  the layer is displayed in the wrong position (way off).

Question:

Why does GeoServer use these transformation parameters and not the official ones, and where does it get them from?

Why don't the official parameters work with GeoServer?